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STATUS  OF  RAY  THFORY  1 WELOPMENT  AT  NAVAL  UNDERSEA 
RESEAT AND  DEVELOPMENT  CENTER 

by 

D.F.  Gordon 

Naval  Undersea  Research  ,na  Development  Center 
San  Dieto,  California,  U.S. 


INTRODUCTION 

The  first  pt  •  of  this  paper  discusses  NUC  (Naval  Undersea  Research 
and  Development  Center)  work  on  the  accuracy  and  validity  of  ray 
theory.  By  comparing  the  range  to  convergence  zones  as  indicated 
by  experiment  and  by  theory,  we  have  found  which  techniques  are 
required  to  make  accurate  ray  computations. 

By  comparing  computations  done  by  ray  theory  and  normal--mode  or 
wave  theory,  we  can  determine  the  limits  to  the  accuracy  of  ray 
theory  at  low  frequencies. 

The  final  portion  of  this  paper  discusses  new  developments  in 
ray  theory. 

CONVERGENCE  ZOSE  V.NGE 

The  range  to  the  leading  edge  of  the  convergence  zone  can  be 
determined  very  accurately  experimentally.  This  is  partly  because 
the  average  travel  time  to  the  leading  edge  of  the  zone  is  very 
insensitive  to  minor  variations  in  the  velocity  profile  and  can 
be  used  to  measure  range  accurately,  and  partly  because  of  the 
rather  abrupt  increase  in  sound  pressure  at  the  zone.  The  range 


at  which  the  propagation  loss  decreases  to  less  than  95  dB  has  been 
used  to  indicate  the  leading  edge  of  the  zone  because  at 
frequencies  of  a  few  kiloHertz  this  loss  is  clearly  distinguishable 
from  bottom  returns.  It  is,  thei  ifore,  interesting  to  compare 
this  »  „;ge  with  the  range  to  the  first  zonal  caustic  of  ray  theory. 

Pedersen  and  Anderson  gave  a  paper  on  this  topic  at  the  28th  Naval 
Symposium  on  Underwater  Acoustics.  Figure  1  is  a  summary  of 
portions  of  that  paper.  The  figure  indicates  average  results  from 
a  number  of  Pacific  locations.  The  jagged  line  indicates  a  possible 
experimentally  observed  convergence  zone  edge.  The  three  vertical 
lines  represent  computed  losses  at  caustics  with  their  characteristic 
shape  and  indicate  the  range  relative  to  the  true  zone.  This  leading 
caustic  is  formed  by  rays  which  travel  downward  from  the  source  and 
upward  to  the  receiver. 

Early  attempts  to  compute  the  range  to  this  caustic,  used  sound 
velocities  computed  from  Kuwahara’s  tables.  Several  trials  gave 
ranges  which  averaged  1.5  kyd  short  of  the  zone,  as  indicated  in 
Fig.  1.  Altering  the  profile  to  simulate  the  effect  of  earth 
curvature  shortened  the  range  an  additional  600  yd  as  shown. 

The  advent  of  Wilson* s  equations  for  the  computation  of  sound 
velocity  increased  the  range  to  computed  caustics  and  the  correction 
for  earth  curvature  became  an  asset  rather  than  a  liability  The 
caustic  line  labelled  "Wilson"  indicates  the  average  relative  position 
of  computed  caustics  from  12  different  locations  in  the  Pacific. 

This  average  position  is  a>out  500  yd  beyond  the  true  convergence 
zone.  However*,  by  applying  a  caustic  correction  taken  from 
Brekhovukikh,  the  difference  between  theory  and  experiment  is  reduced 
by  half. 


The  final  step  in  Pedersen  and  Anderson’s  investigation  was  to  make 
an  adjustment  in  sound  velocities  to  fit  a  portion  of  Wilson’s  data 
which  includes  those  ranges  of  temperature,  salinity,  and  pressure 
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found  in  the  Pacific.  This  gave  sound  velocities  which  were  smaller 
at  shallow  depths  and  larger-  at  greater  depths  than  those  obtained 
from  Wilson’s  equations  by  amounts  up  to  1  ft/s.  The  median 
difference  between  computed  and  experimental  convergence  zone 
ranges  became  zero.  This  is  shown  by  the  diffraction  curve  labelled 
"Wilson  adjusted".  These  results  indicate  in  a  statistical  sense 
that  the  current  methods  for  computing  convergence  zone  ranges  have 
no  significant  bias.  Although  the  average  of  the  differences  is 
zero,  their  scatter  is  not.  Fifty  percent  of  the  differences  were 
less  than  360  yd.  The  largest  difference  was  over  4  kyd. 


RAY  AND  NORMAL-MODE  THEORY 


Figure  2  shows  two  profiles  of  special  form  for  which  both  ray  and 
normal-mode  computations  can  be  made.  Comparisons  will  be  made 
between  computationr  made  by  the  two  theories.  On  the  left  is  an 
Epstein  layer  which  is  a  five-parameter  function  of  hyperbolic  cosines 
and  tangents.  It  has  been  fitted  to  an  Indian  Ocean  velocity  profile. 
The  curve  has  two  vertical  asymptotes,  one  at  I636  and  one  at  1753  yd/s. 
To  simplify  certain  aspects  of  the  problem,  computations  were  done 
without  che  surface,  so  the  profile  is  shovn  extending  above  the 
surface. 


On  the  right  side  of  the  figure  is  a  feur-layer  approximation  to 
an  Atlantic  profile  in  which  the  squared  index  of  refraction  is 
linear  in  each  layer.  Normal-mode  computations  for  this  profile 
at  10  Hz  and  30  Hz  were  published  by  Tolstoy  and  Clay  in  JASA  in  i960. 


Figure  3  compares  propagation  loss  as  computed  in  three  different 
ways  for  the  Epstein  layer.  Since  no  surface  or  bottom  is  included 
in  this  profile  model,  only  energy  trapped  in  the  duct  by  diffraction 
is  seen  at  the  zones.  Only  two  caustics  appear  at  each  zone.  With 
surface  reflection,  throe  additional  caustics  would  appear.  The 
source  and  receiver  are  at  depths  of  33  yd  and  100  yd  and  the 
frequency  is  30  Hz.  The  channel  axis  is  at  1589-yd  depth. 


The  normal-mode  theory  gives  the  most  accurate  solution  for  this 
idealized  duct.  The  inability  of  the  simple  ray  theory  to  compute 
diffraction  effects  is  apparent.  In  the  modified  ray  theory 
Brekhovskikh * s  caustic  correction  has  been  applied  to  each  caustic 
and  the  results,  which  are  Airy  functions,  have  been  added  in  random 
phase.  A  possible  explanation  for  the  difference  between  the  mode 
and  modified  ray  theory  results  is  that  the  caustic  corrections 
were  not  added  in  phase. 

The  next  four  figures  will  compare  ray  and  mode  theory  for  the  four- 
layer  Atlantic  profile.  Mode  theory  will  be  shown  for  10  Hz, 

30  Hz,  and  100  Hz.  Figure  4  is  a  ray  diagram  for  a  source  at  500-yd 
depth  and  the  upper  500  yd  is  shown.  Rays  are  drawn  at  each  1°  in 
source  angle  with  the  rays  that  just  penetrate  into  the  surface 
duct  and  just  graze  the  bottom  included.  The  range  is  to  100  kyd 
and  includes  one  convergence  zone . 

The  leading  caustic  runs  from  54  upward  to  62-kyd  range  before  it 
encounters  the  surface.  A  similar  caustic  is  formed  by  the  rays 
which  start  upward  at  the  source.  Three  additional  caustics  are 
formed  by  the  surface-reflected  rays,  the  last  between  70  kyd  and 
73  kyd. 

Figure  5  shows  propagation  loss  contours  as  computed  by  normal  modes 
for  precisely  the  same  situation  as  was  used  on  the  ray  diagram, 
except  that  an  100  yd  in  depth  is  shown.  The  frequency  is  10  Hz 

The  two  refracted  caustics  and  the  final  surf ace-ref lected  caustic 
from  the  previous  figure  are  shown  by  broken  lines.  Note  that  the 
leading  caustic  from  the  ray  diagram  from  54  kyd  to  57  kyd 
approximately  parallels  the  80  and  90-dB  contours.  Note  also  the 
surf ace -image  effect  which  depresses  the  90-dB  contour  in  the  zone 
deeper  than  50  yd  from  the  surface. 

The  110-dB  contour  appears  to  be  influenced  by  the  surface  duct 
which  has  a  depth  of  153  yd.  However,  judging  from  the  next  figure, 
this  is  not  a  result  of  the  surface  duct  which  is  too  small  at  this 
frequency  to  have  any  large  effect  upon  the  loss. 


Figure  6  is  the  same  as  Fig.  5  but  for  30  Hz.  Here  the  leading 
caustic  lies  partly  within  the  80-dB  contour.  The  second  refracted 
caustic  lies  near  the  string  of  80-dB  contours.  The  90-dB  contour 
comes  within  about  25  yd  of  the  surface  at  this  frequency. 


At  IG-Hz  frequency  the  110-dB  contour  in  the  near  field  extended 
to  28  kyd  range.  Here  it  reaches  only  21  kyd.  It  seems  more 
reasonable  to  attribute  this  difference  to  differences  in  diffraction 
into  the  shadow  zone  than  to  attribute  it  to  propagation  in  the 
surface  duct  which  should  be  stronger  at  the  higher  frequency. 

Figure  7  shows  the  contoured  field  at  100  Hz.  Here  definite  surface 
duct  propagation  is  seen.  This  surface  duct  can  trap  one  mode  at 
100  Hz  so  this  propagation  is  not  unexpected.  The  effect  of  the 
surface  duct  c.,.n  be  seen  in  the  90,  100,  and  110-dB  com  ours  in  and 
following  the  direct  field  and  in  the  100  and  110-dB  contours 
following  the  zone.  The  zone  itself,  as  outlined  by  the  90-dB 
contour,  is  only  slightly  larger  than  the  ray  theory  zone  bounded 
by  the  first  refracted  caustic  and  the  last  surface-reflected  caustic 
between  54  and  73-kyd  range. 

Some  phase  interference  or  Lloyd-mirror  beats  can  be  seen  near  10  kyd. 
They  were  at  somewhat  shorter  range  for  the  lower  frequencies. 

These  figures  have  shown  several  limitations  of  ray  theory  at  low 
frequencies.  Diffraction  from  caustics  and  shadow  zones  is  important 
as  is  the  interaction  between  ducts  such  as  the  S0FAR  duct  and 
surface  duct.  This  interaction  between  ducts  can  remain  important 
at  higher  frequencies.  The  surface  image  or  surface  decoupling  effect 
must  be  considered. 

NEW  TECHNIQUES 
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Three  items  under  current  development  at  NUC  are  generalized 
velocity  functions,  two-dimensional  velocity  variation,  and  numerical 
quadrature . 
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In  March  1968,  Pedersen  published  his  generalized  ray  theory 
in  JASA.  This  theory  uses  depth  as  a  function  of  velocity  to 
represent  the  velocity  profile.  The  function  can  be  a  polynomial, 
a  power  series,  or  a  series  in  non- integral  powers  of  velocity. 
This  makes  it  possible  to  fit  velocity  profiles  directly  with 
polynomials  of  any  required  degree  or  to  use  standard  profile 
forms  by  expanding  velocity  as  a  power  series  in  depth  and  then 
inverting  the  power  series.  By  using  non-integral  powers  of 
velocity,  Pedersen  was  able  to  develop  a  theory  of  the  axial  ray 
published  in  JASA  in  January  1969.  This  ra\  theory  requires  the 
use  of  elliptic  integrals.  However,  new  developments  have 
determined  the  range  and  travel  time  as  a  power  series,  making 
elliptic  integrals  unnecessary.  A  report  by  Pedersen  and  White  on 
this  development  was  given  at  the  recent  International  Acoustic 
Congress  in  Budapest. 


In  another  new  development,  White  and  Keir  at  NUC  have  developed 
a  method  of  determining  ray  fields  with  two-dimensional  velocity 
variations.  This  is  done  by  transformations  on  the  depth  and 
range  axes.  This  technique  gives  theoretical  examples  of  two- 
dimensional  velocity  variation  which  can  approximate  various 
realistic  situations  and  also  can  give  models  to  test  numerical 
ray  tracing  programs. 


In  May  1971,  Mr  Edward  R.  Floyd  of  NUC  published  an  article  in  JASA 
on  ray  tracing  by  Gaussian  quadrature.  This  method  again  allows 
a  polynomial  of  arbitrary  degree  to  be  fitted  to  all  given  velocity 
points  and  thereby  avoids  false  caustics.  It  is  not  yet  clear 
whether  this  method  can  give  sufficient  accuracy  for  computing 
intensities  from  detailed  velocity  profiles.  However,  it  appears 
to  be  well  suited  for  quick  approximations  to  range  and  travel  time, 
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SUMMARY 


The  range  -to  convergence  zones  can  be  accurately  computed  if 
accurate  velocity  profiles  which  are  independent  of  range  are 
known,  and  if  earth  curvature  and  diffraction  from  the  caustic 
are  taken  into  consideration. 

The  technique  of  comparing  ray  and  mode  solutions  for  identical 
velocity  profiles  gives  valuable  information  on  the  validity 
of  ray  theory  for  finite  wavelengths. 

New  work  includes  power  series  expansions  for  a  general  class  of 
velocity  profiles,  velocity-depth  transformations  to  simulate 
two-dimensional  velocity  variation,  and  numerical  quadrature. 


DISCUSSION 

Bartberger  had  also  encountered  convergence  difficulties  using 
Gaussian  quadrature  even  with  25  points.  The  author  felt,  however, 
that  numerical  techniques  were  now  available  which  might  make  the 
method  usable. 

In  reply  to  a  question  concerning  the  continued  use  of  the  random- 
phase  addition  of  modes,  the  author  said  that  certain  results  to  be 
found  in  Brekhovskikh' s  work  now  made  this  unnecessary. 
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A  THEORETICAL  METHOD  FOR  IRE  PREDICTION 
OF  UNDERWATER  EXPLOSION  PULSES  AT  CAUSTICS 


by 

I.M.  Blatstein 
(read  by  R.M.  Barash) 

Naval  Ordnance  Laboratory 
White  Oak,  Silver  Spring,  Maryland,  U.S. 


Our  concern  [Refs.  1  and  2]  is  with  the  effect  of  refraction  on 
the  long  range  propagation  of  underwater  explosion  shock  waves. 

Here,  as  with  acoustic  sources,  ray  tracing  can  be  used  to  predict 
refraction  effects.  From  the  divergence  or  convergence  of  raysr  an 
amplification  factor  can  be  calculated.  This  is  defined  as  the 
square  root  of  the  ratio  of  the  cross  sectional  area  between  rays 
at  a  given  point  assuming  spherical  spreading  to  the  cross  sectional 
area  between  rays  at  the  same  point  when  the  actual  sound  velocity 
profile  is  specified.  We  can  then  multiply  the  pressure  history 
expected  at  a  given  point  if  no  re fraction  occurred  by  the 
appropriate  amplification  factor.  This  then  gives  us  the  pressure 
history  expected  at  that  point  when  refraction  is  accounted  for. 

However,  the  amplification  factor  is  inversely  proportional  to 
the  square  root  of  the  cross  sectional  area  between  adjacent  rays. 

So  as  we  approach  a  caustic,  where  these  rays  cross,  the 
amplification  factor  reaches  infinity,  and  ray  theory  is  invalid. 
Furthermore,  in  the  shadow  zone  adjacent  to  the  caustic,  ray  theory 
predicts  zero  energy  penetration.  This  is  due  to  the  high  frequency 
nature  of  the  ray  ti  eory  approximation.  So  if  we  are  interested 
in  the  pressure  near  a  caustic  or  in  an  adjacent  shadow  zone  a 
method  other  than  ray  theory  must  be  used. 
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In  this  paper,  I  will  describe  such  a  method  for  calculating  shock 
wave  pressure  histories  in  and  near  caustic  regions.  This  method 
involves  the  incorporation  of  various  propagation  effects  into  a 
Fourier  series  representation  of  the  initial  shock  wave  from  an 
underwater  explosion.  I  will  then  describe  comparisons  that  have 
been  made  between  calculated  pressure  histories  and  experimental 
results  from  ocean  and  flooded  quarry  tests.  Figure  1  shows  a 
typical  ray  diagram  for  the  ocean  case  which  we  have  considered. 

Here  the  source  explosion  is  at  a  depth  of  1000  ft.  The  convergence 
zone  caustic  then  occurs  at  20  mi  to  30  mi  from  the  source.  If  the 
source  is  deeper,  in  or  just  below  the  thermocline,  a  thermocline 
related  caustic  occurs.  This  shows  up  at  2  mi  to  5  mi  and  is  due 
to  upward  starting  rays.  The  flooded  quarry  test  that  we  have 
considered  was  intended  to  model  this  thermocline  related  caustic. 

For  both  cases,  comparisons  will  be  shown,  agreement  will  be  praised, 
and  discrepancies  will  be  sullenly  discussed.  Finally,  I  will  talk 
about  the  accuracy  of  ray  theory  near  the  caustic,  a  region  wheve 
it  is  known  to  break  down. 

A  starting  point  for  our  work  was  a  solution  to  the  wave  equation 
that  has  been  done  at  various  times  in  slightly  different  forms  by 
Brekohvskikh  [Ref.  3-],  Tolstoy  [Ref.  4],  and  in  this  case,  Sachs 
and  Silbiger  [Ref.  5]*  [Fig*  2],  The  wave  equation  is  solved  for  a 
harmonic  source  and  an  arbitrary  depth  dependent  sound  velocity 
profile.  First  they  arrive  at  a  ray  solution,  where  each  term  in 
the  sum  corresponds  to  an  arrival  reaching  the  point  of  interest 
after  leaving  the  source  at  a  different  initial  angle.  But  as  with 
any  ray  solution,  this  one  breaks  down  on  the  caustic.  So  they  do 
a  further  approximation  and  arrive  at  an  expression  valid  at  and  near 
a  caustic  for  a  sinusoidal  acoustic  source.  If  we  divide  this 
expression  by  the  pressure  expected  at  that  point  assuming  spherical 
spreading,  we  get  an  amplification  Factor  valid  at  a  caustic  [Fig.  3]* 

Most  of  the  quantities  in  the  amplification  factor  are  constants  of 
the  propagation  path.  These  are  determined  once  the  sound  velocity 


profile,  source  depth,  and  point  of  interest  are  given.  The 
expression  also  contains  the  Airy  function,  a  function  of  both 

3/ 

frequency  through  k'3  and  distance  off  the  caustic  through  h r. 

For  a  given  source  frequency,  the  Airy  function,  and  so  the 
amplification  factor,  falls  off  exponentially  as  we  move  to  the  left 
of  the  caustic  into  the  shadow  zone.  ‘As  we  start  moving  to  the  right 
of  the  caustic,  which  is  located  at  hr  —  0,  the  amplification  factor 
at  first  increases.  Further  to  the  right  we  see  an  oscillating 
function.  This  is  the  result  of  two  arrivals  interfering  in  what  is 
familiar  to  us  in  ray  theory  as  the  double  arrival  region.  From  an 
asymptotic  expression  for  the  Airy  function,  it  can  be  shown  that 
one  arrival  is  approaching  the  caustic  and  has  been  amplified.  The 
other  arrival  is  receding  from  the  caustic.  It  has  been  amplified 
and  also  phase  shifted  by  90°.  While  this  information  about  phase 
shifts  and  amplitudes  is  difficult  to  verify  in  the  re'.ultant  signal 
from  a  sinusoidal  source,  it  is  more  readily  seen  in  the  resultant 
signal  from  a  transient  such  as  a  shock  wave.  We  will  see  this  in 
some  of  the  figures. 

Now  that  we  have  a  frequency  dependent  amplification  factor,  we  need 
to  apply  it  to  a  shock  wave.  Our  next  step  is  to  describe  the  shock 
wave  in  a  frequency  dependent  manner.  This  has  to  be  done  in  such  a 
way  that  we  describe  what  the  shock  wave  would  look  like  at  the  range 
of  the  caustic  if  no  refraction  occurred.  Then  we  can  incorporate 
our  amplification  factor  into  the  pressure  history  in  order  to 
account  for  refraction. 

We  represent  the  initial  shock  wave  as  an  abrupt  rise  to  a  peak 
pressure  followed  by  an  exponential  decay  [Fig.  4]-  We  then  write 
the  Fourier  series  for  the  pulse.  This  way  it  is  expressed  in  a 
frequency  dependent  manner.  In  this  figure,  the  amplitude  of  the 
shock  wave  is  normalized  to  one.  In  general,  we  will  need  both  a 
peak  pressure,  Pi,  and  a  decay  constant,  9,  to  determine  the  pulse 
shape  and  series.  We  would  like  these  parameters  to  be  characteristic 
of  the  range  of  interest  and  also  to  take  into  account  the  finite 
amplitude  effects  that  are  present  in  the  propagation  of  a  shock 
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wave . 


We  do  this  by  using  the  similitude,  or  scaling,  equations  which  were 
discussed  in  the  companion  paper  by  Barash  and  doertner  [Ref,  6] 

(see  Sect.  3  of  these  Proceedings).  Since  these  equations  are  based 
on  cube  root  scaling,  they  properly  account  for  the  finite  amplitude 
effects  on  the  peak  pressure  and  decay  constant.  However,  they 
have  not  been  verified  out  to  the  30  mi  range  at  which  we  need  Pi 
and  3.  So  we  use  them  to  a  range  where  they  are  known  to  be  valid, 
and  the  peak  pressure  is  low  enough  so  that  we  may  ignore  finite 
amplitude  effects  beyond  that  range.  The  point  at  which  we  cease 
to  account  for  finite  amplitude  effects  is  the  range  at  which  the 
peak  pressure  drops  below  5  psi.  From  this  point  to  the  range  of 
the  caustic  we  assume  acoustic  spherical  spreading  to  find  Pi 
and  9  for  our  pulse. 

We  now  have  the  shock  wave  in  a  frequency  dependent  form  with  the 
appropriate  peak  pressure  and  0.  The  next  step  is  to  incorporate 
our  amplification  factor  into  the  expression  [Fig.  5]. 

The  expression  at  the  top  of  the  figure  is  the  Fourier  series  for 
the  pressure  history  with  refraction  added.  Pi  and  0  are  the  pulse 
parameters  we  have  just  discussed.  The  function  represented  by 
the  script  fp  to  the  right  of  the  summation  sign  is  the  refractive 
amplification  factor.  This  along  with  the  Ti/4  added  to  the 
arguments  of  the  sine  and  cosine  take  into  account  the  effect  of 
refraction  on  the  propagating  pulse  near  the  caustic.  This  new 
expression,  if  evaluated  as  it  stands,  would  diverge.  The  physically 
unrealizable  step  discontinuity  in  the  pulse  would  lead  to  an 
infinite  spike  when  refraction  was  added.  However,  we  have  yet  to 
take  into  account  the  attenuation  of  acoustic  pressure  disturbances 
that  becomes  important  for  long  ranges  and  high  frequencies.  Not  only 
is  it  an  important  effect  over  the  long  propagation  path  to  the 
caustic,  but  it  serves  to  force  our  series  to  converge  and  terminate 
at  some  finite  frequency.  The  graph  in  the  centre  of  this  figure 
shows  the  relative  strength  as  a  function  of  wave  number,  k,  of  the 
refraction  factor  f,  and  an  absorption  factor  P.  We  can  see  that 
the  ocean  acts  as  a  filter,  damping  the  high  frequency  components 
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of  the  shock  wave.  Thus  by  adding  absorption,  we  get  the  expression 
at  the  bottom  of  this  figure  for  the  pressure  near  a  caustic.  It  is 
a  Fourier  series  that  now  terminates  at  some  finite  frequency  due 
to  absorption.  This  limitation  shows  up  as  the  summation  now  only 
extends  co  N  rather  than  infinity..  Once  we  find  the  range  to  the 
caustic  for  the  depth  of  interest,  and  determine  the  absorption 
cut-off,  N,  the  expression  for  the  pressure  may  be  thought  of  as 
being  purely  a  function  of  Ar,  the  distance  off  the  caustic.  The 
next  figures  show  pulses  calculated  for  various  Ar's,  along  with 
experimental  pressure  histories  from  the  convergence  zone  experiment 
we  treated  [Fig.  6a].  The  pulse  on  the  right  is  a  typical  experimental 
record.  Those  on  the  left  are  calculated.  For  negative  Ar,  we 
are  in  the  shadow  zone.  Here  wc  see  a  broad,  low  amplitude  pulse. 

This  is  a  result  of  the  action  of  the  Airy  function,  which  is 
monotonically  decreasing  with  frequency  in  this  region.  In  the 
caustic  region  [Fig.  6b],  we  observe  a  single  arrival  with  a  high 
amplitude  and  short  decay  time.  As  we  move  away  from  the  caustic, 
the  peak  pressure  at  first  grows  and  then  starts  to  decrease.  Also 
shown  in  the  figure  is  the  isovelocity  pressure  history.  This  is  what 
the  pressure  history  would  look  like  at  the  range  of  the  caustic  if 
no  refraction  occurred.  For  1 arger  distances  off  the  caustic 
[Fig.  6c],  there  are  two  arrivals,  the  resultant  of  the  two  arrivals 
from  various  frequencies.  The  first  has  just  been  amplified  and 
has  approximately  the  same  decay  constant  as  the  isovelocity  pressure 
history.  The  second  arrival  las  been  amplified  and  further  phase- 
shifted  by  90° .  This  phase  shift  leads  to  a  much  shorter  decay 
constant  as  seen  in  the  figure.  Thus,  for  a  shock  wave,  the  effect 
of  the  phase  shift  shows  up  clearly. 

At  this  point,  I  have  described  the  method  for  calculations  near  a 
caustic  and  shown  qualitatively  what  to  expect.  Next,  I  will  show 
the  comparisons  of  calculated  pressure  histories  and  experimental 
results  that  have  been  made. 

The  first  such  comparison  involves  an  oceanic  experiment  done  by 
the  Naval  Ordnance  ,  >oratory. 


The  purpose  of  this  experiment  was  to  record  pressure  histories  in 
a  convergence  zone.  Figure  1  showed  an  average  sound  velocity 
profile  and  associated  ray  diagram  for  the  time  period  of  the 
experiment.  During  the  experiment,  one  ship  set  off  charges  of 
8  lb  and  900  lb  of  TNT  in  a  region  from  near  the  surface  to 
1000  ft  deep.  Another  ship  some  30  mi  away  used  a  vertical  array 
of  more  than  100  hydrophones  to  record  pressure  histories  throughout 
the  convergence  zone .  When  the  gauge  array  was  placed  so  that  it 
crossed  the  caustic  for  a  given  shot,  pressure  histories  were  then 
obtained,  for  the  shadow  zone,  caus-*  ic  region,  and  double  arrival 
region,  simultaneously. 

We  selected  five  of  these  shoiti  tc  test  the  method  just  described. 

We  will  show  representative  pressure  histories  from  four  of  these 
shots . 

In  order  to  do  the  calculations,  we  must  know  the  depth  of  the 
gauge  string  with  respect  to  the  caustic.  For  these  shots,  this 
was  determined  from  an  analysis  of  the  experimental  pressure 
histories.  In  Fig. 7,  the  position  of  the  gauge  string  for  each 
shot  is  indicated  on  an  enlarged  view  of  the  upper  region  of  the 
convergence  zone.  We  notice  that  there  is  a  reflected  branch  of 
the  caustic  resulting  from  surface  reflection,  as  well  as  the 
direct  branch.  Gauges  at  various  positions  will  record  different 
pressure  histories  due  to  their  proximity  to  both  branches  of 
the  caustic.  For  the  first  two  shots  considered,  due  to  lack  of 
complete  data,  we  will  only  consider  the  contribution  from  the 
direct  caustic. 

We  will  first  consider  shot  151  [Fig.  8],  an  8  lb  shot  where  the 
gauge  string  was  relatively  deep.  For  shot  151,  we  first  specify 
the  gauge  of  interest,  for  example,  the  gauge  near  the  top  of  the 
gauge  string.  Then  we  calculate  the  parameters  needed  for  the 
amplification  factor  on  the  caustic  at  that  depth  using  the  given 
sound  velocity  profile  and  source  depth.  Then  from  this  figure  we 
determine  A r,  t;-.v  horizontal  distance  from  the  gauge  to  the 
caustic,  and  the  final  parameter  needed. 
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Figure  8  contains  this  calculated  pressure  history,  and  the 
experimental  record  from  this  gauge  in  the  shadow  zone  in  the  upper 
left-hand  corner.  Next  to  it  is  a  comparison  in  the  caustic  region, 
and  at  the  bottom  of  the  slide  are  two  comparisons  in  the  double 
arrival  region.  In  all  cases  the  experimental  record  is  a  solid 
line,  while  the  calculated  record  is  a  dashed  line.  In  all  three 
regions,  the  calculated  pressure  histories  adequately  match  the 
experimental  records.  In  the  upper  right-hand  corner  we  have  also 
plotted  the  isovelocity  pressure  history.  This  again  is  what  the 
pressure  history  would  look  like  at* all  four  gauges  at  the  range 
of  the  caustic  if  no  refraction  occutred. 


Also  of  interest  is  the  good  agreement  for  relatively  long  times 
aft >r  the  peak  pressure.  The  caustic  solution  used  in  our  pressure 
history  retains  some  of  the  high  frequency  limitations  of  ray  theory. 
This  would  suggest  that  for  each  arrival  we  could  only  make  valid 
predictions  near  the  peak  pressure  whei  j  high  frequencies  predominate, 
Yet  these  calculations  yield  good  results  beyond  the  point  where  the 
pressure  history  drops  below  zero  and  flattens  out.  This  indicates 
&  reasonably  good  description  of  the  low  frequency  content  of  the 
pulse.  This  is  probably  due  to  the  relatively  small  sound  velocity 
gradients  in  the  ocean  which  make  the  caustic  solution  valid  for 
relatively  low  frequencies  on  the  order  of  llO  Hz  or  less.  This  low 
frequency  validity  should  extend  to  ray  theory  as  well,  despite  the 
often  made  remark  that  ray  theory  is  valid  for  high  frequencies  only. 


The  next  shot  we  consider  is  shot  82  [Fig.  9],  an  890  lb  explosion 
at  1000  ft.  The  same  general  agreement  between  experimental  records 
'.nd  calculated  pressure  histories  is  evident.  In  the  records  closest 
to  the  caustic,  [Figs.  9b&9c],  gauge  case  ringing  is  very  severe 
and  interferes  somewhat  with  the  comparisons. 


Now  we  will  consider  the  two  shots  where  both  the  direct 
and  reflected  arrivals  were  recorded.  A  gauge  at  the  bottom  of  the 
gauge  string  for  shot  120  is  in  the  double  arrival  region  of  the  direct 
branch  caustic  and  in  the  shadow  zone  of  the  reflected  branch  caustic. 


I 


Hence,  we  would  expect  arrivals  due  to  both  branches.  From  ray 
geometry  we  can  calculate  the  arrival  time  difference,  or  time 
delay,  between  these  arrivals  from  different  branches.  We  further 
assume  that  the  surface  acts  as  a  perfect  reflector  for  the 
reflected  arrivals,  causing  only  a  phase  reversal.  We  then  combine 
the  resulting  pulses  with  the  appropriate  time  delay  to  find 
the  pressure  history. 


Figure  10  shows  comparisons  of  these  calculated  pressure  histories 
and  experimental  records  from  the  top,  middle,  and  bottom  of  the 
gauge  string.  In  the  bottom  record,  the  pressure  history  starts 
with  a  direct  double  arrival,  which  is  then  followed  by  a  negative 
shadow  zone  pulse  from  the  reflected  caustic.  Again  the  general 
agreement  of  peak  pressure  and  wave  forms  is  good. 


As  the  last  of  the  oceanic  comparisons,  we  calculated  pressure 
histories  for  shot  119,  an  8  lb  shot  [Fig.  11].  This  differed 
from  the  previous  shot  in  that  it  is  a  smaller  charge  weight,  and  so 
we  are  dealing  with  pulses  with  a  much  shorter  decay  time .  This  means 
there  is  proportionately  more  high  frequency  energy.  By  treating 
records  with  reflected  arrivals  in  them,  we  are  able  to  test  the 
method  for  considerable  distances  off  the  caustic,  in  Fig,  11,  the 
reflected  arrival  is  approximately  16GG  ft  1 orizontally  into  the 
shadow  zone,  and  the  direct  double  arrival  is  approximately  2400  ft 
horizontally  from  the  direct  caustic.  So  with  one  record  we  gain 
information  about  a  region  4000  ft  wide . 


In  all  four  comparisons,  the  calculated  pressures  are  in  reasonable 
agreement  with  the  experimental  values  from  the  oceanic  test.  Just  as 
important,  the  entire  waveforms  are  in  good  agreement.  So  for  the 
oceanic  convergence  zone  case,  the  method  is  a  reasonable  approximation 
for  the  various  phenomena  involved  in  propagation  to  a  convergence 
zone . 

We  also  tried  matching  pressure  records  from  a  test  conducted  by 
Woods  Hole  * n  a  flooded  quarry  [Fig,  12].  As  I  have  said,  this  type 
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of  test  models  the  thermocline  related  caustic  that  may  occur  at 
2  mi  to  5  mi.  For  the  quarry  case  the  solution  to  the  wave 
equation  is  not  clearly  valid  for  the  frequency  domain  of  our  pulse. 
Not  only  are  the  sound  velocity  gradients  1000  times  larger  in 
the  quarry  than  in  the  ocean,  but  the  caustic  is  much  closer  to 
the  turning  points  of  the  rays.  This  could  tend  to  restrict  the 
validity  to  the  higher  frequencies  in  the  pulse. 


p 3 


Also,  we  have  to  modify  the  pressure  expression  since  absorption 
is  no  longer  the  high  frequency  cut-off  mechanism  as  it  was  in  the 
ocean.  Here  the  gauge  response  restricts  high  frequencies  more 
than  absorption,  so  the  gauge  response  as  a  function  of  frequency 
is  used  in  the  pressure  expression  instead  of  absorption. 

Pressure  histories  were,  then  calculated  near  the  caustic  and  in 
the  double  arrival  region  for  a  56  lb  charge  detonated  at  50  ft 
[Fig.  13].  The  dashed  lines  are  the  calculated  pulses,  the  solid 
lines  are  the  experimental  pressure  records.  While  the  exact 
agreement  of  the  peak  pressures  is  no  doubt  fortuitous,  good 
agreement  between  experimental  and  calculated  peak  pressures  has  been 
obtained  for  similar  quarry  data  during  the  course  of  this  analysis. 
However,  it  is  also  apparent  from  this  figure  that  the  decay  from 
the  peak  pressure  is  much  too  steep. 

If  the  solution  to  the  wave  equation  is  indeed  valid  only  for  high 
f:  equencies  for  this  case,  it  may  explain  why  our  calculations  only 
appear  reasonable  near  the  peak  pressure  where  high  frequencies 
predominate.  Also,  it  is  possible  that  a  propagation  mode  other  than 
the  pure  refracted  one  contributes  to  the  pressure  an  a  thermocline 
related  caustic.  For  example,  a  lateral  wave  type  of  propagation 
could  be  occurring  at  the  interface  of  the  almost- isovelocity  upper 
layer  in  the  flooded  quarry.  This  would  be  energy  at  relatively  low 
frequencies  which  would  affect  most  the  decay  from  the  peak.  Or  it 
is  possible  that  the  surface  boundary  is  influencing  the  propagation 
of  energy  in  a  way  not  accounted  for  by  ordinary  ray  theory.  If  this 
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were  the  case,  a  modified  ray  theory  such  as  that  described  by 
Murphy  and  Davis  in  these  Proceedings  might  be  necessary  for  this 
type  of  situation.  Whatever  the  reason,  the  calculation  of  the 
decay  of  the  pulse  for  the  quarry  case  remains  the  aspect  of  this 
work  that  requires  the  most  improvement. 

Finally,  I  would  like  to  discuss  the  validity  of  ray  theory  near 
the  caustic.  Qualitatively  we  know  that  ray  theory  predicts  an 
increasing  intensity  as  one  approaches  the  caustic  from  the  double 
arrival  region.  This  intensity,  which  is  the  same  for  all 
frequencies,  reaches  infinity  on  the  caustic,  a  line  of  zero  width. 
The  caustic  solution  that  we  are  using  yields  a  different  picture. 

For  this  solution,  the  caustic  is  a  region  of  finite  thickness  in 
between  the  shadow  zone  and  double  arrival  region.  The  lower  the 
frequency  of  the  source,  or  the  larger  the  charge  for  the  transient 
case,  the  wider  the  caustic  region  is.  Inside  the  caustic  region, 
the  solution  yields  one  arrival  that  slowly  increases  in  amplitude 
as  one  moves  away  from  the  shadow  zone  boundary  at  &r  =  0.  For  large 
distances  off  the  caustic,  the  expression  yields  two  arrivals  in  the 
pressure  history.  This  is  the  double  arrival  region  as  expected 
from  ray  theory.  However,  since  the  caustic  solution  was  derived 
to  be  valid  near  the  caustic,  it  is  not  obvious  that  it  should  yield 
reasonable  results  far  from  the  caustic  in  the  double  arrival  region. 
So  really  there  are  two  questions  to  be  answered:  How  different  are 
ray  theory  and  the  caustic  solution  in.  the  caustic  region  where  ray 
theory  is  supposed  to  be  breaking  down?  and  How  do  ray  theory  and 
the  caustic  solution  compare  in  the  double  arrival  region,  where  the 
caustic  solution  was  not  originally  intended  for  use? 

Referring  to  Fig.  14,  for  reasons  of  simplicity,  calculations  were 
done  for  a  harmonic  source  rather  than  the  pulse  sources  we  have 
been  treating.  For  the  figure  we  have  assumed  a  100  Hz  harmonic 
source  and  the  convergence  zone  sound  velocity  profile  previously 
shown.  In  the  caustic  region,  we  have  combined  the  two  ray  arrivals 
in  order  to  compare  the  resultant  to  the  caustic  solution  we  have 


been  using.  As  expected,  as  we  approach  the  caustic,  the  two 
solutions  diverge  as  the  ray  solution  increases  rapidly.  It  should 
be  remembered  that  this  figure  holds  true  for  the  particular 
convergence  zone  profile  treated.  We  have  examined  others  where 
the  caustic  region  was  narrower,  and  some  where  the  caustic  region 
was  up  to  three  times  as  wide  as  bhia  one  is.  Also  shown  in  the 
figure  is  the  double  arrival  region.  Here  we  can  resolve  two 
arrivals  in  the  caustic  solution.  It  turns  out  that  each  one 
has  the  same  amplitude,  which  falls  off  as  the  fourth  root  of  Ar. 

We  have  plotted  this  amplitude  along  with  those  of  the  two  ray 
arrivals.  If  we  combined  the  caustic  solution  arrivals  at  the 
caustic  region  boundary,  the  resultant  would  match  the  caustic 
region  solution  at  the  caustic  region  boundary.  However,  what  is 
most  striking  is  the  agreement  between  ray  theory  and  the  caustic 
solution  far  into  the  double  arrival  region.  This  explains  why  our 
calculations  of  pulses  in  the  double  arrival  region  were  in  good 
agreement  with  the  experimental  pulses.  While  in  the  double  arrival 
region,  one  would  tend  to  trust  ray  theory  more,  this  figure 
indicates  that  it  is  not  a  bad  approximation  to  continue  using  the 
caustic  solution  that  w«s  used  in  the  shadow  zone  and  caustic,  region. 

In  order  to  put  all  of  these  results  in  their  proper  perspective, 
we  must  keep  in  mind  the  assumptions  and  approximations  made . 

We  assumed  that  the  unrefracted  shock  wave  could  be  approximated  by 
an  abrupt  rise  followed  by  an  exponential  decay.  We  then  wrote  the 
Fourier  series  for  this  pulse.  In  the  oceanic  case,  we  assumed  that 
a  superposition  of  finite  amplitude  effects,  spherical  spreading, 
and  absorption  on  this  series  would  adequately  describe  long  distance 
propagation  in  the  absence  of  refraction.  Then  we  added  refraction 
effects  to  find  the  pressure  near  a  caustic.  In  the  quarry,  the 
gauge  response  was  substituted  for  absorption.  Our  ability  to  apply 
these  effects  separately  is  no  doubt  due  to  the  phenomena  being 
treated.  Finite  amplitude  effects  occur  relatively  close  in  where 
pressures  are  high.  Absorption  for  high  frequencies  takes  effect 
only  over  long  distances.  While  refraction  effects  occur  throughout 
the  path,  the  major  effect  is  the  significant  amplification  at  the 
convergence  zone.  So  while  these  effects  no  doubt  interfere,  their 
major  influence  is  individual. 
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The  main  thrust  of  this  work  has  been  to  describe  the  effect  of 
refraction  on  the  propagation  of  an  underwater  explosion  shock  wave. 
Particularly,  we  are  interested  in  calculations  near  the  caustic 
where  the  refractive  amplification  is  greatest.  We  feel  that  in 
this  method  we  have  a  reasonable  means  of  calculating  just  such 
pressure  histories  for  the  oceanic  convergence  zone.  What  I  also 
hope  has  been  demonstrated  by  this  paper,  and  the  companion  one 
by  Barash  and  Goertner  [Ref.  6],  is  that  an  underwater  explosion 
can  be  more  than  just  a  source  of  energy  at  various  frequencies 
for  acousticians  interested  in  harmonic  sources.  These  transient 
sources  can  also  lead  to  analyses  and  information  about  ray  theory 
validity  that  would  be  otherwise  difficult  to  obtain. 
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DISCUSSION 


During  the  discussion  the  points  were  made  that  in  a  sense  one 
was  "saved  by  absorption",  that  taking  more  terms  in  the  Fourier 
series  expansion  would  not  help,  and  that  the  use  of  a  mix  of 
linear  and  non-linear  properties,  although  expedient,  may  not  be 
valid  at  large  distances. 


SOUND  VELOCITY  PROFILE  AND  RAY 
DIAGRAM  FOR  A  CONVERGENCE  ZONE  CAUSTIC 
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SAMPLE  VELOCITY  PROFILE  AND  RAY  DIAGRAM  INDICATING 
WAVE  EQUATION  SOLUTION  PARMEiERS  (AFTER  SOBER) 
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THE  REGION  OF  DEFINITION  FOR  THE  EXPONENTIAL  DECAY 
REPRESENTED  BY  A  FOURIER  SERIES 
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PRESSURE  HISTORY  WITH  REFRACTION 
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EXAMPLES  OF  PRESSURE-TIME  HISTORIES  FOR  VARIOUS 
HORIZONTAL  DISTANCES  FROM  A  CONVERGENCE 
ZONE  CAUSTIC 
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EXAMPLES  OF  PRESSURE-TIME  HISTORIES  FOR  VARIOUS 
HORIZONTAL  DISTANCES  FROM  A  CONVERGENCE 
ZONE  CAUSTIC 
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EXAMPLES  OF  PRESSURE-TIME  HISTORIES  FOR  VARIOUS 
HORIZONTAL  DISTANCES  FROM  A  CONVERGENCE 
ZONE  CAUSTIC 
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SHOT  120,  W=870  LB.  COMPARISON  OF  DATA  AND  CALCULATED  PRESSURE  HISTORIES 
(A)  2g  =  44FT.,  Ar2=  B04FT.;  Ar,  =  -91FT.;  (B)  2,  =  121FT.,  Af2=  1219FT. 

Ar,  =  -693FT.;  (C)  Zg  =  202FT.,  Ar2=  1S81FT.,  Ar,  =  -1565FT. 


SHOT  111,  W4  LB.  COMPARISON  OF  DATA  AND  CALCULATED  PRESSURE  HISTORIES. 
(A)  2,-2  FT.,  Af2M2«  FT..Af,-3IO  FT.;  (B)  2j=M  FT..  Ar2=1082  FT. 
Ar,=-254  FT.;  [C)  Zj=202  FT.,  Ar2-2365  FT.,  Ar,=-1565  FT. 
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VELOCITY  PROFILE  AND  RAY  DIAGRAM  FROM  BLACKINGTON 
FARM  QUARRY  FOR  DAY  OF  WHOI  CASE  #45 


FIG.  12 
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INTENSITY  AT  CAUSTICS 


C.W.  Spofford 

Bell  Telephone  Laboratories 
Whippany,  New  Jersey,  U.S. 


The  classical  divergence  expressions  predict  an  infinite  intensity 
at  caustics.  This  infinite  intensity  is  a  fundamental  limitation 
of  geometrical  acoustics  (ray  theory).  The  intensity  at  a  caustic 
is  of  primary  concern  since  the  caustic  is  assumed  tc  be  the  point 
of  highest  intensity  within  a  convergence  zone.  Several  approaches 
to  estimating  the  field  at  a  caustic  have  been  made  using  stationary 
phase  and  more  refined  asymptotic  methods. 


The  recent  uniform  asymptotics  work  of  Ludwig  and  Kravtsov  may  be 
used  to  obtain  both  the  intensity  at  the  caustic  and  the  transition 
from  the  caustic  to  the  point  where  ray  theory  is  applicable. 
Whereas  in  the  context  of  uniform  asymptotics,  coherent  ray  theory^ 
is  the  zeroth  order  solution  (in  wavenumber)  to  the  wave  equation, 
Ludwig  shows  that  the  first  order  solution  is  given  (excluding  a 
phase  factor)  by 

'Kr)~'n1/8  j(Ax  +  A2)  x1/*  Ai  (-x)  +  i  (As  -  A\  )  x1/4  Ai  »  (-x)  }  [Eq.  1] 


where 


x  =  [|  uj(ts  -  Tx)]! 


[Eq.  2] 


Coherent  ray  theory  estimates  intensity  by  adding  all  amplitudes  on  a  phased 
basis. 


ami  Ai,  As,  Ti ,  and  T2  are  the  amplitudes  and  travel  times  of 
the  two  rays  intersecting  the  point  of  interest  (R)  on  the 
illuminated  side  of  the  caustic.  tt)  is  the  angular  frequency, 

Ai  and  Aif  are  the  Airy  function  and  its  first  derivative. 

This  result  may  be  used  to  study  the  validity  of  coherent  ray 
theory  near  a  caustic.  As  one  moves  away  from  the  caustic  (on  the 
illuminated  side)  ^  goes  uniformly,  with  wavenumber,  to  the  result 
obtained  from  coherent  ray  theory.  In  fact,  if  'll  is  compared  with 
the  coherent  sum  of  the  paths  (adding  a  -tf/2  phase  shift  to  Ray  2 
after  tangency  to  the  caustic),  one  finds  good  agreement  up  to 
the  last  point  (first  point  in  range)  of  constructive  interference. 
After  this  point,  the  coherent  sum  tends  to  infinity,  whereas  the 
uniform  asymptotics  result  experiences  an  exponential  decay  into 
the  shadow  zone. 

The  intensity  at  the  last  peak  is  3*5  dB  greater  than  that  at  the 
caustic  and  the  peak  is  frequently  displaced  from  the  caustic  by 
many  wavelengths.  Also  the  position  of  this  peak  corresponds  to 
the  last  point  of  constructive  interference  as  given  by  ray  theory, 
and  the  amplitude  is  given  to  within  0.5  dB  by  the  in-phase  sum  of 
the  geometric  amplitudes. 


The  field  at  the  caustic  proper  is  given  by 

V2  [uu(Tm  R"  -  T'^"’)]1/6 


'li~22/aTTVfe  Ai  (0) 


°R  cos  es 


cs  sin  9rr 


(R")s/J 


[Eq.  3] 


where  R”,  Rw ,  T"  and  Tm  are  the  second  and  third  derivatives  of 
range  and  travel  time  with  respect  to  9g  at  fixed  yR,  and  ®R 
is  the  angle  at  the  receiver.  These  derivatives  may  be  evaluated 
directly  from  the  ray  trace.  In  fact,  for  vertically  stratified 
media,  for  which  Snell’s  law  holds,  'll  may  be  written  in  terms  of  R". 
The  results  of  Fig.  1,  however,  suggest  that  it  may  not  be  necessary 
to  evaluate  Eq.  3  for  the  intensity  at  a  caustic  but  merely  calculate 
the  coherent  sum  of  the  geometrical-acoustics  amplitudes  at  the  last 
point  of  constructive  interference . 
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ACOUSTIC  FIELD  AT  CAUSTICS 


UNIFORM  ASYMPTOTICS 
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COHERENT  RAY  THEORY  ( - )  VS  UNIFORM  ASYMPTOTICS  (-) 


SPECIAL  FORMULATION  OF  MODIFIED  RAY  ANALYSIS 
FOR  MACHINE  COMPUTATION 


E . L .  Murphy 

SACLANT  ASW  Research  Centre 
La  Spezia,  Italy 


J.A.  Davis 

Woods  Hole  Oceanographic  Institution 
Woods  Hole,  Mass-,  U.S. 


When  the  reflection  coefficient  for  waves  incident  on  a  refracting 
or  bounded  region  has  a  phase  h( 0o)  that  io  a  function  of  the  angle 
of  incidence  60  (measured  at  some  reference  level  Z  =  Z0),  then 
in  consequence,  individual  rays,  or  "beams"  will  be  displaced  from 
the  location  predicted  by  Snell's  law.  [See,  for  example, 

Ref.l,  Ch.  I),  By  h(0o)  we  mean  phase  changes  other  than  those 
corresponding  to  the  geometrical-optics  phase  integral  along  the 
ray  path. 

This  phenomenon  of  beam  displacement  is  the  basis  for  modified 
ray  analysis  considered  in  this  paper.  The  development  of  the 
theory  is  given  in  Refs.  2  to  5  (the  nomenclature  is  somewhat 
altered  for  convenience  from  the  symbols  used  in  these  references) . 


The  problem  in  obtaining  quantitative  results  lies  in  the 
determination  of  the  reflection  coefficient 


R(90)  =  |P.|  e 


in  ( 0O) 


[Eq.  1] 


for  realistic  sound  speed  profiles.  In  particular,  the  modification 
of  ray  analysis  are  especially  significant  for  rays  with  turning- 
points,  or  vertices,  near  boundaries,  or  near  a  local  maximum  in 
sound  speed. 

To  determine  R( 0O)  we  use  an  extension  of  the  geometrical-optics 
approximation  in  a  form  that  can  be  made  valid  throughout  regions 
near  a  vertex  on  a  ray.  This  procedure  does  leave  some  freedom 
in  the  choice  of  the  sound  speed  profile  c(z)  other  than 
those  for  which  exact  solutions  are  known. 

The  method  is  such  that  the  phase,  arg  R (  Q 0)  =  k(80)>  and  the 
amplitude  |r(0o)|,  and  many  features  of  the  modified  ray  analysis 
can  be  described  entirely  in  terms  of  a  parameter  E(0O)  determined 
by  the  integral, 


n  p  Z  2  (  8  0  )  . . . — 

E(80)=i--k0  Jn  (z)  -  sin  80dz 

zr(e0) 


[Eq.  2] 


where  n(z)  =c0/c(z)  is  the  index  of  refraction  referred  to  a 
reference  sound  speed  c0  =  c(z0).  The  angle  80  is  measured  with 
respect  to  the  vertical,  or  z-direction,  at  a  point  on  the  ray  at 
reference  depth  z  =  z0.  The  symbol  k0  is  a  reference  wave  number, 
k0=2nf/c0,  where  f  is  the  frequency.  (We  are  considering  a 
formalism  for  a  harmonic  source,  but  we  can  extend  the  application 
of  the  results  to  other  time  functions  by  the  use  of  Fourier 
analysis).  The  limits  of  the  integral  are  zeros,  zi  and  za,  of 
the  integrand  for  a  given  value  of  8,,.  The  reason  for  the 
appearance  of  two  zeros,  zi  and  z2,  when  a  ray  in  reality  has  only 
one  real  vertex,  say  at  z =  zi ,  is  indicated  in  the  sketch  in  Fig.  1 
It  turns  out  that  modified  ray  results  are  most  significant  when 
c(z),  or  n(z),  is  such  that  there  are  at  least  two  zeroes  of  the 
integrand  in  Eq.  2.  Therefore  in  this  paper  we  consider  profiles 
with  a  local  maximum  in  sound  speed,  or  where  a  boundary  at  z  =  z^, 


as  sketched  above,  gives  rise  to  an  "image”  of  the  vertex  at  z  -  zx 
We  confine  our  discussion  to  the  three  kinds  of  problems: 

1)  Zeros  near  a  local  maximum  in  sound  speed 
(unbounded  problem) . 

2)  Zeros  near  a  pressure  release  boundary. 

3)  Zeros  near  a  rigid  boundary. 

The  results  for  the  phase  k(E)  of  the  "plane  wave"  reflection 
coefficients  for  these  three  problems  are  as  follows  [see  Ref.  4]. 


x(E)  =  -‘/.(E)  -  ~  +  e  arctan  e 


[Eq.  3] 


e  =  0  unbounded 

=  -1  pressure  release  boundary 
=  +1  rigid  boundary 


m 


^  *  ’ 


where  the  spatial  function  x(E)  is  given  by  the  expression 


X(E)  - 1  -  |  ^  -l|J-  +  Im  Mr  (-1~2-1-E)  ♦ 


[Eq.  4] 


The  symbol  refers  to  the  natural  logarithm  of  the  absolute 

value  of  E.  The  last  term  is  the  imaginary  part  of  the  logarithm 
of  the  gamma  function  [see  Ref.  6,  Ch.  6]. 

In  fig.  2,  these  phases  are  plotted  as  a  function  of  the  parameter  E. 
Large  positive  values  of  this  ray  Darameter  E  correspond,  for  the 
bounded  problems,  to  rays  that  hit  the  surface,  and  the  phases 
approach  the  l-‘mits  we  usually  associate  with  such  boundaries  when 
the  bounded  medium  is  homogeneous  (-tt  for  pressure  release,  0  for 
rigid  boundary) .  On  the  other  hand,  large  negative  E-values 
correspond  to  rays  with  vertices  far  from  the  boundaries,  and  all 
cur1  as  for  phase  approach  the  valn«  -tt/2,  the  phase  change  for  plane 
waves  reflected  in  a  smooth  monotonic  profile. 


The  curves  are  continuous  and  not  the  simple  "patched-up"  step 
function  that  would  result  if  we  assume  n  =  -tt  for  all  rays 
hitting  the  surface,  or  -n/2  for  all  rays  with  vertices. 

The  amplitude  of  the  reflection  coefficient  is,  of  course,  unity 
for  the  bounded  problems, while  for  the  unbounded  problem  it  has 
the  following  relatively  simple  form, 


|r(e)  1  - 


(1  +  e"c') 


ttEv  l/a 


[Eq.  5] 


The  formula  for  the  displacement  AR,  of  a  ray  that  has  transited 
a  "two-turning-point"  region  is  determined  from  the  derivative 
of  k  with  respect  to  sin  0O  [see  Ref.  2], 


tsR  -  ~  rr~ 


dn 


dE  dn 


k0  d  sin  80  k0  d  sin  0O  dE 

The  derivative  dx/dE  can  be  expressed  as  follows? 


[Eq.  6] 


j 

! 


_  d*  J2  sech  TIE 
dE  dE  e  4  4 


where  dy/dE  can  be  written  in  the  fo 


rm 


SX  =  .ifcM  +  I 

dE  2  2 


j.  „  ,  /I  +  iE  % 

2  Re  V  ^  2  > 


[Eq.  7] 


[Eq.  8] 


The  last  term  involves  the  real  part  of  the  digamma  function  \Sf 
[see  Ref,  6,  p.  258]. 

Since  h(E),  dx/dE,  and  |R{E,)|  are  functions  of  E  only,  then 
the  preparation  of  computer  programs  for  the  contribution  of  these 
functions  to  modified  ray  analysis  is  relatively  simple.  We  only 
need  to  program  their  evaluation  as  functions  of  E.  When  the  actual 
sound  speed  profile  is  selected  then  E ( 8 0)  and  the  derivative 
dE/d(sin  60)  needed  for  AR  can  be  '*valuated  separately. 


£ 

d 

4 


In  our  current  programs,  which  have  been  written  for  the  Woods  Hole 
Oceanographic  Institution’s  SIGMA  7  computer,  and  Hewlett-Packard 
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computer,  a  composite  linear/parabolic  profile  for  l/c2 ,  or 
equivalently,  for  n2(z),  is  used.  The  function  n2 (z)  and  its 
first  derivative  are  made  continuous  at  some  "matching"  distance  P, 
as  sketche,d  in  Fig.  3* 

In  the  sketch  at  the  bottom  of  Fig.  3,  for  a  more  or  less  realistic 
profile,  the  regions  in  which  "two-turning-point"  phenomena  may 
occur  are  shaded.  Perhaps  these  programs  may  eventually  be  useful 
as  simple  subroutines  added  to  already  existing  ray  programs. 

When  some  "flag"  tells  a  ray  it  has  transited  su^h  a  region,  the 
subroutine  supplies  a  displacement.  The  ordinary  ray  program  in 
its  usual  way  then  continues  the  displaced  ray  until  it  again  reaches 
a  "diffraction"  region.  If  the  diffraction  region  is  a  region  of 
local  sound  speed  maximum  the  subroutine  also  supplies  a  change  in 
amplitude  due  to  3.eakage  and  splitting  of  rays. 

The  functions  x(E)  and  dy/dE  are  the  only  terms  in  x(E)  and 
dx/dE  requiring  special  procedures  for  machine  computation. 


With  series  expansions  for  the  gamma  and  diagamma  functions  [see 
Ref.  6,  Ch.  6]  inserted  into  Eqs,  4  and  8,  |x|  and  dx/dE  can  be 

put  into  the  form 


i  X I  =• 


and 


g=-|to2|E!-^+  lEf 


m: 


=1,3,5..  m(m2  +  |e(2) 


[Eq.  9] 


[Eq.  10] 


where  y  Is  Euler’s  constant, 

y  0.57721  _ _ 


The  lo -rarithmie  terms  in  Eqg  ,  9  and  10  lead,  to  computational 
difficulties  since  they  become  infinite  for  E =  0,  that  is  for 
the  grazing  ray  in  the  bounded  case,  or  the  split  ray  in  the 
unbounded  case,  however,  these  infinities  rather  than  being 
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troublesome  are  one  of  the  most  physically  significant  contributions 
of  the  modified  ray  analysis.  The  following  illustration  shows  why. 


Consider  a  purely  parabolic  profile  for  n  (z),  for  example. 


n2  (z)  =  a  +  bz2  . 


[Eq.  lx] 


For  a  source  and  receiver  at  depths  z0  below  the  profile  extremum, 


the  ordinary  ray  theory  result  for  the  range  R  of  the  ray 

r  ay 


emitted  at  angle  80  and  connecting  source  and  receiver,  is  given 
by  the  relation 


1/3 


2  sin  0 


o  Oo-  zl]  +  Zo 


ray 


i/Ti 


Z1 


[Eq.  12] 


where  zx  is  the  depth  of  the  ray  vertex  below  the  extremum  of  the 
n2-profile.  Note  the  logarithmic  infinity  for  the  grazing  or 
split-ray  (z1=0).  However,  for  the  profile  of  Eq.  11,  the  integral 
for  E  can  be  evaluated  to  give, 


E  =  ~k0  J b  z2 
and,  in  turn. 


[Eq.  13] 


dE 


2  sin  8, 


k  -  d  sin  8, 


[Eq.  14] 


In  consequence,  it  follows  that  the  logarithmic  term  in  the 
displacement  is  of  the  form 


2  sin  80  i/a  i/4 

&R  =  - — -  &i(k0  b  '  z\  )  +  other  terms  . 

J  b 


[Eq.  15] 


Therefore,  in  the  modified  range  defined  by  the  relation 


^MQD  Rray  ’’’  *'A 


[Eq.  16] 


the  logarithmic  singularities  cancel  leaving  a  finite,  wavelength- 
dependent  range.  For  machine  computation,  therefore,  we  remove 
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the  logarithmic  singularities  from  the  displacement  AR  and  Rray 
before  programming. 

This  feature  of  the  modified  ray  theory  that  would  arise  for  all 
smooth  profiles  (continuous  derivatives  in  the  region  of  the  profile 
extremum),  leads  to  the  result  that  there  is  a  split-beam  shadow 
at  finite  range  when  modifications  are  included.  Without  the 
modifications  every  point  is  illuminated  by  a  ray. 

On  the  other  hand,  for  a  bilinear  profile  (discontinuous  derivatives) 
ordinary  ray  theory  does  give  a  split-b«am  shadow  at  finite  range, 
that  is  to  say,  there  is  no  singularity  in  Rray«  For  "non-smooth" 
profiles  (discontinuous  derivatives)  we  can  also  develop  a  modified 
ray  analysis  [see  Ref.  4]«  It  turns  out  that  in  this  case  there  is 
no  singularity  in  the  displacement  AR.  In  fact,  as  we  shall  see, 
when  the  modified  ray  results  are  compared  for  a  bilinear  profile  and 
a  smooth  profile  where  the  smoothing  compared  with  the  bilinear  case 
is  only  over  a  depth  increment  of  a  few  wavelengths,  the  results  are 
approximately  the  same.  The  modifications  introduced  by  modified 
ray  theory  have  removed  some  of  the  pathological  sensitivity  of 
ordinary  ray  analysis  to  small  changes  in  the  sound  speed  profile. 

For  the  presentation  of  the  results  of  the  computer  program,  we 
currently  use  the  format  sketched  in  Fig.  4.  The  cycle  range  R( &) 
for  source  and  receiver  at  the  same  depth  is  plotted  as  a  function 
of  the  grazing  angle  6  measured  at  the  source  point  on  a  ray. 

In  Fig.  5  the  solid  curves  represent  examples  for  the  three  kinds 
of  problems  (1)  unbounded,  (2)  pressure  release  boundary,  and 
(3)  rigid  boundary.  In  the  small  figure  in  the  upper  left-hand  corner, 
the  shape  of  the  R  vs  6  curve  based  on  unmodified  ray  analysis  is 
sketched;  of  course  it  is  the  same  for  all  three  problems,  but  different 
for  the  bilinear  and  composite  profiles.  Note  the  singularity  in 
the  ordinary  ray  theory  results  for  a  smooth  profile.  Note  also 
that  for  the  ordinary  ray  theory  result  fo*’  the  bilinear  profile  the 
extremum  in  range  would  not  be  considered  a  caustic,  in  fact  the 
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derivatives  of  the  R  vs  5  curve  are  discon+  '  nuous  for  this 
extremum  (grazing  ray) . 

The  following  conclusions  can  be  drawn  from  the  results  presented 
in  Fig.  5: 

1)  The  modifications  have  removed  the  singularity  in  the 
R  vs  6  curves  for  smooth  sound  speed  profiles.  The  resulting 

vs  6  curves  are  smooth  so  that  the  extremum  does  represent 
a  caustic. 

2)  The  modifications  also  lead  to  smooth  Rj^qd  vs  &  curves 
for  the  bilinear  profile,  so  again  the  extremum  actually  is  a 
caustic. 

3)  The  derivatives  of  the  vs  &  curves  can  be  found 

and  used  as  parameters  in  quantitative  descriptions  (Airy  integrals) 
for  the  field  in  the  neighbourhood  of  this  caustic.  Pedersen  [Ref.  7] 
has  shown  how  falsa  caustics  can  arise  in  ordinary  ray  analysis  for 
linear-segmented  profiles.  It  is  rather  interesting  to  see  that 

i 

with  modified  ray  analysis,  the  modifications  not  only  may  act  to 
remove  some  of  this  pathology,  but  also  may  result  in  real  caustics 
at  places  not  recognized  as  caustics  in  the  ordinary  ray  analysis 
for  bilJnear,  or  bounded  linear  n  (z)  profiles. 

4)  The  caustics  are  displaced  in  range  compared  with  the 
location  we  might  anticipate  on  the  basis  of  ordinary  ray  analysis. 

For  rays  with  vertices  near  a  pressure  release  boundary,  displacements 
are  toward  the  source;  for  a  rigid  boundary,  displacements  are  mostly 
away  from  the  source. 

5)  The  modifications  have  made  the  results  different  for 
all  three  kinds  of  problems,  even  for  rays  whose  vertex  would  lie 
below  the  region  where  the  models  differ. 

6)  For  the  bilinear  profile,  and  for  a  bilinear  profile 
"smoothed"  over  a  depth  increment  of  "half-width"  P  =  2  wavelengths, 
the  results  are  approximately  the  same . 
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In  the  lower  part  or  Fig.  5,  the  values  of  the  magnitude  of  the 
reflection  coefficient  for  the  unbounded  problem  are  plotted  to 
give  quantitative  results  for  the  splitting  of  the  rays  into 
reflected  and  transmitted  branches.  The  diffraction  regions  can  be 
sources  of  ’’leakage"  rays  and  diffraction  reflected  rays  (from  rays 
that  according  to  ordinary  ray  theory  would  be  totally  transmitted) . 

In  some  of  the  examples  we  have  analysed,  it  appears  that  for 
frequencies  of  the  order  of  a  kiloHertz  the  displacements  can  be 
hundreds  of  feet;  while  at  lower  frequencies,  of  100  Hz  or  so, 
displacements  amount  to  thousands  of  feet. 

Finally,  we  would  like  to  point  out  that  whereas  here  we  have 
considered  a  harmonic  source  represented  as  an  integral  over  ray 
parameter  0O,  it  is  also  possible  to  consider  signals  of  different 
time  behaviour,  Fourier-analysed  so  that  the  integral  is  over  the 
variable  ui.  In  that  case,  we  can  also  show  that  there  are  time 
delays  or  advances  At,  in  arrival  times,  that  are  related  to 
derivatives  of  the  phase  with  respect  to  ui.  It  can  be  shown 
that  these  time  displacements  can  also  be  written  in  terms  of  the 
derivative  dx/dE,  as  follows, 


.  _  dx  _  dx  dE 
^  du)  dE  cTU) 


[Eq.  17] 


E 

U)  dE 


Recently,  a  number  of  laboratories  have  been  using  machine  programs 
for  modal  analysis  where  a  sufficient  number  of  modes  (perhaps  a 
thousand,  if  necessary)  are  summed  to  describe  the  sound  field  for 
comparison  with  ray  computations.  It  should  be  borne  in  mind,  when 
making  such  comparisons,  that  the  modal  analysis  is  a  more  complete 
and  exact  analysis  for  the  field,  and  if  there  is  any  validity  for 
the  displacements  or  modifications  we  have  described  in  this  paper, 
the  modal  results  should  already  include  such  effects. 
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DISCUSSION 


In  answer  to  questions  the  first  author  said  that,  broadly  speaking, 
modifications  to  ray  analysis  were  needed  if  the  ray  has  a  vertex 
within  4  or  5  acoustic  wavelengths  of  a  surface}  and  that  the 
linear/parabolic  model  with  "half -width"  p  was  chosen  for  purposes 
of  illustration. 

The  first  author  affirmed  that  amplitudes  could  be  obtained 
straight-forwardly  from  modified  ray  analysis.  In  reply  to  a 
suggestion  that  this  was  equivalent  to  including  higher  order 
terms  in  an  inverse-wavenumber  expansion,  the  first  author  said 
that  this  may  be  so,  but  that  the  individual  terms  were  not  the 
same  as  in  the  traditional  expansion. 

Some  discussion  at  this  session,  and  in  later  private  sessions, 
revealed  that  the  techniques  of  modified  ray  analysis  might  be 
usefully  applied  to  remove  some  of  the  pathological  aspects  of 
liscontinuities  when  velocity  profiles  are  approximated  by 
constant-gradient  segments. 
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THE  EFFECT  OF  GRAVITY-FORCED  OSCILLATIONS  AT  THE  BASE 
OF  THE  DUCT  ON  ITS  EFFECTIVE  DEPTH 
AS  A  CHANNEL  FOR  ACOUSTIC  RAYS 

by 

I.  Roebuck 
(read  by  G.  Murdoch) 

Admiralty  Underwater  Weapons  Establishment 
Portland,  Dorset,  U.K. 


This  talk  must  begin  with  a  justification,  since  at  first  sigh 
(and  possibly  even  after  further  glances)  it  is  difficult  to  s 
what  place  a  problem  of  classical  applied  mathematics  (which  was 
originally  chosen  as  a  companion  to  the  author’s  earlier  paper  on 
below  duct  sound  levels  due  to  scattering  from  the  sea  surface)  can 
rightfully  claim  to  occupy  in  this  gathering  of  ray  tracers  and 
computing  enthusiasts.  This  is  especially  so  when  one  realises  that, 
as  in  the  earlier  paper,  given  here  six  months  ago,  the  object  is  to 
formulate  the  problem  in  such  a  way  as  to  obtain  the  maximum  of 
predictive  results  with  the  minimum  of  computation. 

Nevertheless,  this  account  of  the  effect  of  gravity-forced  internal 
waves  at  the  base  of  an  isothermal  layer  on  its  efficiency  as  a  duct 
for  acoustic  energy  is  opposite  to  this  particular  convention;  it 
has  the  laudable  aims  of  saving  effort  for  the  sceptical 
perfectionists,  and  saving  face  for  the  lazy  intuitionists,  by 
demonstrating  that  it  is  not  necessary  to  develop  three-dimensional 
ray-tracing  programs,  nor  to  solve  numerically  Ihe  wave  equation  for 
a  duct  with  one  sinusoidal  boundary,  in  order  tc  take  account  of 
the  horizontal  stratification  of  sound  velocity  introduced  by 
internal  waves. 
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This  is  a  very  worthwhile  simplification,  but  to  achieve  it  it  is 
necessary  to  have  some  background  as  to  the  generation  and  acoustic 
effects  of  internal  waves;  we  start,  as  always,  with  the  classical 
picture  of  acoustic  propagation  in  the  upper  ocean,  with  the 
isothermal  surface  duct,  and  associated  underlying  shadow  zone  at 
long  ranges  [Fig.  1],  This  is  the  ideal,  on  which  computations 
of  field  intensities  can  be  made  with  relative  ease,  as  the 
stratification  purely  in  horizontal  planes  makes  a  two-dimensional 
ray  analysis  valid. 

Unfortunately,  in  the  real  world,  things  aren't  that  simple.  Even 
if  we  retain  the  basic  framework  of  ray  theory,  there  are  three 
major  perturbing  effects  to  be  taken  into  account  when  considering 
the  idealised  situation  as  a  predictive  model  for  in-duct 
propagation.  The  first  of  these  is  diffraction  which,  as  shov/n 
schematically  in  Fig.  2,  can  distribute  energy  into  zones  forbi 
by  simple  ray  acoustics.  This  is  obviously  important  as  far  as 
energy  levels  in  the  shadow  zone  are  concerned;  there  is  a  school 
of  thought  that  claims  that,  as  far  as  in-duct  propagation  is 
concerned,  the  effect  is  negligible.  The  argument  runs  that  the 
effect  of  diffraction  is  merely  to  alter  the  effective  depth  of 
the  duct,  and  that  the  "below  layer"  field  is  an  integral  part  of 
the  trapped  modes.  Even  so  it  is  clear  that  there  is  a  net  energy 
leakage  by  this  mechanism. 

The  second  complication  to  be  considered  is  scattering  from  the 
rough  sea  surface  of  energy  which  is  seemingly  entrapped  within  the 
duct,  so  that  it  is  deflected  into  the  thermocline  region  and 
escapes.  This  is  shown  schematically  in  Fig,  3j  much  vork  ha  been 
carried  out  on  this,  both  numerical  and  analytic  (the  papers  by 
van  Ness,  Schweitzer  and  the  present  author,  to  name  but  three), 
and  it  is  fair  to  say  that  the  effect  of  this  mechanism  is  now 
quantifiably  determined,  or  determinable,  for  most  'typical'  duct 
conditions . 
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Finally,  among  these  mechanisms,  we  come  to  anomalous  refraction 
due  to  internal  waves,  the  least  well  documented  and  most 
speculative  of  the  three.  Schematically,  the  way  in  which  acoustic 
energy  which  has  been  trapped  in  the  duct  can  be  guided  out  of  it 
by  the  action  of  internal  waves  is  shown  in  Fig.  4$  the  dominant 
effect  is  the  change  from  upward  to  downward  curvature  on  crossing 
the  internal-waves  profile.  It  was  on  the  basis  of  this  type  of 
diagram  that  Schulkin  made  his  widely  accepted  estimate  of  the 
effect  of  internal  waves,  reducing  the  < "'fective  duct  depth  by 
the  rms  height  of  the  typical  internal  wavej  the  argument  running 
that  any  ray  which  vertexes  at  a  greater  distance  than  this  from 
the  surface  will  in  time  intersect  such  an  internal  wave  and  be 
refracted  out  of  the  duct. 


Unfortunately,  this  conclusion  is  invalid,  as  this  diagram  is 
totally  misleading  —  an  inevitable  consequence  of  the  distortion 
of  vertical  and  horizontal  scales  to  get  the  figure  on  to  a 
conventional  slide.  What  in  fact  happens,  because  the  curvature 
of  the  acoustic  rays  is  so  small  and  they  are  being  propagated 
almost  horizontally,  is  ohat  the  ray  path,  even  taking  account  of 
refractive  differences,  occupies  several  periods  of  the  internal 
wave  in  any  transition  from  in-duct  to  below-duct  propagation.  The 
true  schematic  is  more  like  the  one  shown  in  Fig.  5  —  again  noting 
that  this  is  grossly  distorted  —  the  true  grazing  angle  to  the 
internal  waves  is  less  than  1°;  this,  however,  at  least  indicates 
tnat  a  ray  may  penetrate  into  the  internal-wave  region  and  even 
so  re-emerge  into  the  surface  duct. 


Now  we  are  getting  to  the  core  of  the  problem  —  if  some  rays  can 
be  refracted  back  into  the  surface  duct  while  others  are  lost  from 
it,  how  do  we  calculate  what  is  the  effect  of  the  internal  waves 
in  quantitative  terms?  To  do  this  obviously  requires  a  more 
detailed  knowledge  of  the  mechanism  of  generation  and  propagation  of 
internal  waves  and  so  I  must  ask  you  to  lay  down  your  acoustics 
and  follow  into  the  uncharted  depths  of  oceanography.  The  forcing 
mechanism  for  internal  waves  is  gravity  (salinity  can  also  cause 


a 

3S 


I 
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them,  but  for  the  purpose  of  this  paper  it  will  be  neglected), 
and  buoyancy  forces  are  dominant. 

We  start  with  a  water  mass  which  was  in  equilibrium  at  a  depth  z 
below  the  surface,  under  a  pressure  p,  its  density  at  that  time 
being  P,  and  which  by  some  mischance  has  been  adiabatically 
displaced  a  small  distance  6z  from  this  position.  The  general 
definition  of  the  word  'adiabatic'  in  this  context  is  beyond  the 
scope  of  this  paper;  in  this  context  it  means  that  if  the  same 
supernatural  agency  that  caused  the  change  in  the  first  place 
decides  to  put  the  water  back  in  its  original  place,  its  density 
and  pressure  will  also  revert  to  their  initial  values.  Even  more 
confusingly,  but  more  importantly,  it  means  that  although  sea 
water  is  a  viscous  non-Newtonian  fluid  we  can  treat  it  as  though 
it  were  a  perfect  gas  and  obeys  the  gas  laws,  Ptt  PT  and  Ptt  p\ 
Anyway,  in  its  new  position,  the  density  of  this  element  is 


whereas  the  local  equilibrium  density  is  P(z+  6z)  .  Thus  this 
element  has  a  density  deficiency  of 


P(z)  6z 


,1  dp  l  dP, 
(P  d z  "  v? 


relative  to  its  surroundings,  and  experiences  a  gravity- fed 
restoring  force.  The  equation  of  motion  is 


w 


I 

I 


That  has  introduced  yet  another  variable,  6,  so  I  had  better 
explain  what  that  is.  It  is  the  temperature  the  element  would 
end  up  at  it'  it  adiabatically  was  transferred  to  a  position  where 
it  was  subject  to  a  standard  pressure  of  1  atmosphere,  and  it  is 
called  the  potential  temperature. 

N  is  a  measure  of  the  speed  of  reaction  of  the  ocean  to  a 
perturbation  —  in  other  words  of  its  stability.  The  bigger  N2 
the  more  stable  the  ocean.  Figure  6  shows  the  actual  distribution 
of  N2  in  the  upper  ocean  in  one  typical  case,  together  with  the 
normal  idealisation  of  it  used  by  theorist **.  Normal  in  this  case 
means 

a.  ti.ct  it  is  the  one  normally  used; 

b.  that  it  is  a  normal  distribution; 

c.  that  any  ocean  to  which  it  is  a  valid  approximation 
is  completely  abnormal. 

So  now  we  know  what  happens  when  a  fluid  element  is  displaced 
vertically;  an  internal  wave,  though,  has  a  horizontal  particle 
velocity  component,  so  the  full  equation  of  motion  must  be  used. 

We  can  simplify  them  by  applying  Boussinesq' s  Approximation,  which 
says  that,  for  slow  enough  motions  (which  internal  waves  are)  we 
can  treat  the  fluid  as  incompressible  except  that  buoyancy  must 
be  taken  into  account.  The  equations  are  the  top  four  in  Fig.  7, 
which  I  do  not  intend  to  go  into  any  detail  over,  except  to  say 
that  you  can  eliminate  all  the  other  dependent  variables  and  come 
up  with  this  general  equation  for  the  vertical  component  of 
velocity,  w 


d3  /  a2  ,  a2 


(_2_  +  JLo  w  +  n2  ^  z )  -  o 


which  holds  no  matter  what  kind  of  internal  motion  we  are  discussing. 
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However,  we  are  after  a  propagating  internal  wave,  so  we  take  the 
most  general  disturbance  which  is  sinusoidally  periodic  in  time 
and  travelling  in  the  horizontal  direction,  writing  [Fig.  8] 

w  =  W(z)  expfi(kx~nb)| 
and  find  that  W  must  obey  the  equation 


dfW  +  N 5  -  ns 


k2W  =  0 


dz‘ 


n 


which  means  immediately  that  the  disturbance  decays  everywhere 
except  inside  the  region  where  its  time  frequency  is  less  than  the 
local  Brunt-  aisala  frequency.  So,  from  Fig.  8  that  shows  the 
distribution  of  N  with  depth,  we  can  say  that  internal  waves, 
except  he  very  slow  ones,  are  confined  to  the  region  of  the 
thermocline.  With  that  fact  established,  we  can  go  further  and 
derive  a  dispersion  relation  between  the  time  frequency  n  and 
the  horizontal  space  frequency  k. 


£  -  kg  _ 

P  1  +  coth  kD 


(kD  »  1) 


In  this  D  is  the  depth  at  which  the  thermocline  is  situated  and 
•bP  is  the  change  in  density  across  it.  Now  in  this,  an  increase  in 
k  causes  an  increase  in  n,  but  n  cannot  be  greater 
than  the  B  r un  t - V a i s a 1 a  frequency,  so  we  can  find  a  maximum  value  of 
k  for  each  value  of  N  .  That  defines  the  minimum  wavelength 

for  a  wave  of  this  frequency;  the  wave's  amplitude  is  limited  by 
the  depth  of  the  region  in  which  the  Brunt-Vaisal’a  frequency  is 
large  enough  to  support  it,  so  we  have  a  relation  between 
amplitude  and  wavelength. 


More  usefully,  we  can  plot  from  Fig,  9,  which  you  can  read  as  either 
the  maximum  slope  a  wave  of  given  amplitude  can  have,  <or  the  maximum 
amplitude  for  a  wave  of  given  slope  on  the  mean  thermocline  plane 
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These  will  ;e  the  ones  which  have  the  greatest  effect  on  acoustic 
propagation,  and  Fig.  10  shows  the  model  adopted  for  calculating 
this  effect.  Sinusoids  make  the  calculation  too  difficult,  so 
they  have  been  replaced  by  truncated  prisms  of  the  same  wavelength 
and  maximum  slope,  cut  off  so  that  the  area  of  ’’intrusion"  is  the 
same  as  that  under  the  sine  wave.  The  assumed  velocity  profiles 
are  the  same  as  the  ones  used  for  the  undisturbed  situation 
with  the  rays  in  the  region  of  intrusion  bei.'g  straight  lines,  as 
the  perturbation  occurs  without  change  of  sound  velocity.  To 
calculate  the  ray  path  we  just  use  Snell's  Law  so  the  calculation 
is  straightforward  but  tedious. 

Because  it  is  so  tedious,  all  that  the  author  has  considered  are 
those  rays  which  would  be  horizontal  at  the  base  of  the  undisturbed 
duct.  He  has  not  finished  even  these,  btt  says  that  no  matter  what 
the  internal  wave,  more  than  90$  of  the  time  the  ray  ends  up  back 
in  the  duct.  So  it  seems  that  everyone  can  go  away  happy  in  the 
knowledge  that  their  previous  neglect  of  internal  waves,  whether 
due  to  ignorance  or  indolence,  was  and  still  is  justified. 
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APPLICATION  OF  THE  RIESZ  POTENTIAL  TO  THE  CAUCHY  PROBLEM 
FOR  WAVE  PROPAGATION  IN  AN  INHOMOGENEOUS  MEDIUM 

t>y 

I .A.  Lopes 

Naval  Undersea  Research  and  Development  Center 
San  Diego,  California,  U.S. 


ABSTRACT 

The  method  of  Riesz  [Ref.  1]  for  the  solution  of  hyperbolic 
partial  differential  equations  is  applied  to  the  Cauchy  problem 
for  the  wave  equation.  It  is  shown  that  the  first  term  in  the 
Riesz  potential  function,  which  is  represented  in  series  form, 
yields  the  geometrical  acoustics  solution  when  applied  to  the 
problem  of  radiation  from  a  point  source. 


THE  WAVE  EQUATION  AND  ITS  RIEMANNIAN  GEOMETRY 

We  seek  a  solution  to  the  partial  differential  equation 

Lu  =  ^-utt-y2u  =  f 


[Eq.  1] 


where  f  is  a  function  of  (t,r)  tR  «  E3  .  For  simplicity  we 
shall  assume  vanishing  initial  conditions,  u(0,r)  =  u^(0,r) =  0. 
The  local  sound  speed  c  is  assumed  to  be  a  function  of 
r  =  (x,y,z)  only. 

Construction  of  the  Riesz  potential  for  the  operator  L  rests 
on  the  Riemannian  geometry  associated  with  the  operator.  The 
semi-Riemannian  metric  is  given  by  the  differential  form 


ds2  =  c2  dt2  -  (dr)J 


[Eq.  2] 
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A  displacement  for  which  ds2  >0  is  called  timelike;  one  for 
;hich  ds2  <0  is  called  spacelike.  Let  P:(to,  rQ)  and 
Q<(t1,r1)  be  two  points  of  space  time.  A  geodesic  joining  P 
and  Q  is  a  curve  y  :  {  [*•  (o ) ,  r  (a )  ]  ;  0<o^o0|  such  that 

[t  (0) ,  r(0]  =  (t0,r0),  rtf0o),  ?(c0)  ]  =  (t1,r1),  and 


s(P,Q)  “  J  ds 
Y 


[Eq.  3] 


is  an  extremum  with  respect  to  all  curves  joining  P  to  Q. 
Let 


(dS)2 

ldaJ 


.2  /dtx2  /dr^ 
'  Ma  Mo' 


=  w 


Then  the  geodesic  curves  are  extremum  curves  of  the  integral 
°o  1 


J  °  W*  do 


f  Eq.  41 


The  geodesic  curves  then  satisfy  the  Euler-Lagrange  equations 

i?  <!£>  -  If1  -  0  r^q.  si 


d^  (*r>  -  v 


w 


2  =  0 


[Eq.  6] 


where 


rO  - 


do 


v  =  ^ 
do 


[Eq.  7] 


[Eq.  8] 


If  we  now  choose  w  to  be  constant  along  a  geodesic  (o  j-s  then 
a  linear  function  of  s)  the  equations  for  the  geodesics 
become 


do 


(c2v°)  = 


vo2c  ac  = 


dt 


0 


=  -(v0)3  c7c 


[Eq.  9] 


[ Eq . 10 j 
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4E2g\H^f?*P 


t 


Let  w  =  1  -  p3 .  Then  p  =  1  is  the  equation  of  the  characteristic 
cone  through  P .  We  have 


/dsv2  _ 

(3j)  =o 


2  _2 

2  yO  _  y 


[Eq.  11] 


If  Oip£l  ds2  >0  and  the  geodesics  are  timelike.  Since  w 
is  constant  along  a  geodesic 


[Eq.  12] 


Using  the  notation  introduced  by  Hadamard  we  set  F(p,Q)  =  s3  (P,Q). 
The  region  |Q:  r(p,Q)  >  0,  t a  tx  \  is  called  the  retrograde 
conoid  with  vertex  P.  We  see  from  Eq.  12  that  it  is  also 
defined  by  0  a  p  <  .  t  £  tx  .  The  region  |Q:  r(P>Q)  ^  0,  t  >  tx  \ 
is  called  the  direct  conoid  with  vertex  P. 


In  the  theory  of  the  Riesz  potential  a  particular  coordinate 
system  for  the  conoid  with  vertex  P  plays  an  important  role. 

This  is  the  Riemannian  coordinate  system  with  coordinates  defined 

by 

X1  =  v^Ojo  ,  i  =  0,  1,  2,  3  [Eq.  13] 

In  this  coordinate  system  the  geodesics  emanating  from  the  point 
P  appear  as  straight  lines.  It  is  shown  by  Riesz  [Ref.  1]  that 

Ar  =  8  +  2ci  [Eq.  14] 

where  A  represents  the  second  order  Beltrami  operator,  or  space- 
time  Laplacian,  and  a  is  the  determinant  of  the  metric  tensor, 
expressed  in  Riemannian  coordinates.  We  h^ve 

a  =  J2  c2  [Eq.  15] 

ft  x  y  z  \ 

where  J  =  D\xo  xi  xq)  >  the  Jacobian  transformation  from 
the  original  coordinate  system  to  the  Riemannian  coordinate 
system.  Expressed  in  our  original  coordinates,  t,  x,  y,  z, 
we  have 
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ar  = 


?-  (i  r  ) 

c  '  c  1  t '  t 


•  (c  7  r) 


—  LF  —  V  0n  c  *  vr 


[Eq.  16] 


=  Lr  +  2a 


d  in  c 
da 


since  ■—  =  0 


Hence,  from  Eqs.  14  and  16, 


LF  =  8  +  2o 


den  J 
do 


f  Eq*  17] 


From  Eq.  9  we  see  that  c2v°  is  constant  along  a  geodesic. 
Letting  the  constant  be  -  c0  =  -c(P)  , 


dt 

da 


[Eq.  18] 


dr 

da 


dv 

do 


~  v 


V  c  = 


1  /  0  \ 
2  V  (~) 


TEq.  19] 

TEq.  20] 


Moreover,  from  Eq.  11, 


‘V(0)]3 


f*  Eq.  21] 


THE  RIESZ  POTENTIAL 

Following  Duff  [Ref.  2]  we  denote  by  D*5  the  interior  of  the 
retrograde  conoid  with  vertex  P.  Let  S  be  the  intersection 
of  D^  with  the  initial  manifold,  t  =  0,  and  let  D*5  be  the 
part  of  the  conoid  cut-off  by  the  initial  manifold,  i.e.  the 
intersection  of  with  the  half-space  t>  0.  For  twice 
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Av~e*-#**.*  «V*=>  T-7*-^'\^»r»*-*Kjpf5>  ,v>?£S\,  P^^TOn^E 


differentiable  functions  u,  v,  defined  on  Dp  we  have  by 

s 

Green's  theorem 

J*  (u  L  v  -  v  L  u)  dt  dx  dy  dz 

DP  [Eq.  22] 

j  “4"  (uv^  “  vut^nt  +  (VXJU  “  UVV)  *  n  dS  1 

p  c 

s 


“X 

sue 


where  CP  is  the  part  of  the  characteristic  cone  (T  =  0) 
s 

cut  off  by  the  initial  plane  t  =  0,  and  (n,  ,n)  is  the 

c 

exterior  normal  to  the  boundary. 

The  Riesz  potential  Va(p,Q)  is  defined  as  a  function  of 
points  P,Q,  in  the  Riemannian  space  of  the  wave  equation 
and  the  complex  variable  a  .  It  satisfies  the  relations 

LVa+2  (P,Q)  =  Va(P,Q)  [Eq.  23] 


and 


lim  P  f(Q)  Va(P,Q)  dt  dx  dy  dz  =  f(P)  ,  [Eq.  24] 

a  -»  0  J 

Dp 

s 


for  any  continuous  function  f.  Ya  is  expressed  in  the  form 

«  sa  +  2k“  \(P,2) 

Va(P’Q)  -  So  . H (a^T -  * 


fEq.  25] 


where  s  =  s(P,Q)  is  the  geodesic  distance  between  P  and  Q. 

The  functions  Y^PjQ)  are  to  be  determined  from  the  conditions 
23  and  24,  while 

H  (a ,  k )  =  tt  2a  +  k_1  r(£a)  f(£a  +  k-l)  [Eq.  26] 

For  sufficiently  large  Rea,  Va  is  an  analytic  function  of  a 

and  vanishes  on  Cp  .  Thus,  for  those  values  of  a,  and  for 

s 

functions  u  satisfying  the  vanishing  initial,  conditions  of  the 
Cauchy  problem  [Eq.  1],  we  have  from  Eq.  22 
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f  u  L  V-a  +  2  dt  dx  dy  dz 
rJ 

DP 


-  r 


yd  +  2  ^u  ^  ^  ^y  ^ 


u 

DP 

s 


[Eq.  27] 


Equation  27  remains  valid  for  all  values  a  to  which  analytic 
continuation  is  possible.  From  Eqs.  23  and  24,  then,  letting 
a  tend  to  zero,  we  obtain  for  a  solution  to  Eq.  1 


u(P)  =  lira  I1  Va  +  2  (P,Q)  f(Q)  dt  dx  dy  dz  .  [Eq.  28] 
a  -»  0  J 


dp 

s 


Equation  28  provides  a  representation  of  the  solution  to  the 
Cauchy  problem.  Determination  of  the  coefficients  Vk(P,Q) 


remains. 


From  the  definition  given  by  Eq.  26  and  the  prooert.ies  of  the 
Gamma  function  we  obtain  the  relations. 


H(p,k)  =  2(i  p  +k  -  2)  H(0,  k  -  1) 


H (P  +  2 ,  k)  =  2p  (g-p+k-l)  H(p,k) 


[Eq.  29] 
f  Eq .  30] 


Now 


L(rp  t )  =  rp  l*  +  *Lrp  1 4p  rp_1  o 


[Eq.  31] 


[Here  we  have  used  Eq.4.5.19  of  Duff  (Ref.  2).]  Using  Eq.  17 
this  becomes 


L(r^)  =  rPL*  +  4prp_1  jo +  (p  +  i  +  -|o^i)  <t«  |.  [Eq.  37- 


do 


32] 


Operating  with  L  on  Vp  +  2(P,Q),  given  by  Eq.  25,  and 


using  Eq.  32,  we  have 


LV 


a  +  2 


£0  H(a  +  2,  k)  i4(2  +  k_1)  r7+K  *  [Eq.  33] 


a 


+  k— 4 


[(|  +  k  +  i  a  )  vk  +  a  -j-|j  +  +  k  "  1  LVk  | 


Using  Eq.  29  this  becomes 


,Ct  +  2 


"  |0  ai^Tinm  I  (a  +  2  >■ +  0^1)  vk 


av  a  .  u  0 

+  2°t£  +  LVk-ll  ^ 


[Eq.  34] 


where  we  have  introduced  V^PjQ)  =  D.  Now  choose  V^P^Q) 
so  that  for  k=0,l,2, 


(2k  +  a  *ML  )  Vk  +  2o  +  L  Vk  _  x  =  0  [Eq.  35] 


Then,  using  Eqs.  29  and  30,  we  have 


LV*  +  2  =  20  inaTITT^lT 


J  +  k-  2 


£0  H“(a  ,kT  Vk  r 


|  +  k-2 


=  V 


[Eq.  36] 


Thus  requiring  Vk  to  satisfy  Eq.  35  results  in  Va  satisfying 
Eq.  23.  We  may  rewrite  Eq.  35  in  the  form 


£(okj*vk>  - 


-iok_1  LVk_1 


[Eq.  37] 


Choosing  V0 (0)  =  c0 ,  so  that  Va  will  satisfy  Eq.  24  [Ref.  1] , 


we  have 


v0(p,g: 


[Eq.  38] 


since  J  *=  1  at  P.  For  ki?l  we  set  V.  (0)  =  0.  Then 


Vk  = 


h  o~k  / c14""1  J*  LVk-i 


[Eq.  39] 


.  ■•  ■.?■;  ■■■  k  ■  "*  > "  ■ .  f  ■!..<  >  .  M%sm 
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APPROXIMATE  SOLUTION  TO  THE  CAUCHY  PROBLEM 


a  ,  2 

Let  us  now  replace  Vu  ^(P,fi)  by  tiie  first  term  in  the  series 
[Eq.  25]  in  computing  the  integral  [Eq.  28],  Define 


V  r!'1 

Vp>  “  J  f<2)  sk+ToT  dtdxiiydz 


[Eq.  40] 


We  consider  ol^n^u^(p)  to  be  a  first  approximation  to  the 

solution  to  the  Cauchy  problem.  In  order  to  carry  out  the 

integration  indicated  in  Eq.  40  we  introduce  a  coordinate  system 

for  D*5  based  on  the  geodesics  defined  by  Eqs.  18,  19,  and  20. 
s 

With  v  =  (v1  ,  v2  ,  v3 )  we  set 


v1  (0) 

v3  (0) 

v3  (0) 


p  COS  6  cos  cp 
p  cos  8  sin  cp 
p  sin  8 


[Eq.  41] 


Then  Eq.  21  is  satisfied.  The  geodesic  equations  provide  a 
correspondence  between  points  (t,x,y,z)  6.  and  (a,p,cp56)» 
This  correspondence  will  not  be  one-to-one  in  general.  If  multi¬ 
paths  occur  a  point  (t,x,y,z)  may  correspond  to  many  points 
(a,p><‘p,0).  However,  except  at  exceptional  points  (t,x,y,z), 
each  of  the  points  (o,p,cp,9)  will  have  a  neighbourhood  that 
is  in  one-to-one  correspondence  with  a  neighbourhood  of 
(t3x,y,z).  At  the  exceptional  points,  called  caustic  or  focal 
points,  this  local  one-to-one  property  will  not  hold.  As  a 
consequence  the  Jacobian  determinant  D(^'x^p  will  vanish  at 
such  points. 


The  geodesic  Eqs.  18,  19  and  20  uniquely  define  (t,x,y,z) 
as  functions  of  (a>p>cp}S).  Transforming  coordinates  in 
Eq.  40  we  have 


2tt  "/  2  1  ,J  t  —  -  1 

J  dcp  j*  d0  J  dp  J  °  d0  f(t,x,y,z)  VQ  [o2  (1  -  p?  )]2  D(^e>) 


[Eq.  42] 


where  0-tQ  =ot  ( p » cp  © )  •*-s  the  value  of  o  for  which  the  geodesic 

curve  reaches  the  initial  manifold.  This  will  happen  at  a  finite 
value  since  we  assume  c  is  bounded.  It  is  well  known  that  if 
F(p,€)  is  a  continuous  function  for  0  <  p  <  1,  £  i  0,  that 


r 


€  J  F(p,€)  (l-p)6"1  dp  -  F(1,0) 


[Eq.  43] 


•a 


-I 


m 


$ 
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We  will  write,  for  continous  g(p), 
Jg(p)  6(1  -  p)  dp  =  g(l) 


[Eq.  44] 


Equation  44  is  the  definition  of  the  generalized  function  s(l-P). 
Since  H(a  +  2,0)  ~4  tt/  a  as  a  -♦  0  we  can  apply  Eq.  43  to  Eq.  42 


to  obtain 


u(p) 


=  lim  u  (?) 
a-*  0  ' 


[Eq.  45] 


2n  n/2  1  ct 

^  j  J  «  j  dp  do  fvo0-a  6(1- p)K^^) 
0  H  o  © 


Note  that  in  the  latter  equation  we  have  assumed  f  V_  d(^  x  ^  to 

n  °  o  p  cp  8 

be  a  continuous  function  of  p  in  a  neighbourhood  of  p  =  1  . 

The  solution  to  the  Cauchy  problem  is  represented  approximately 
b-  .  45.  Multipaths  cause  no  problem  in  this  representation, 
e  they  are  sorted  out  by  the  (o>p>cp>G)  coordinates. 


Lett  inf 


represent  a  point  source  we  set 


i (f  ,x,y,z)  =  S(t)  &(x,y,z) 


[Eq.  46] 


wnera  6(x,y,z)  is  the  3~dimensional  delta  function  and  S(t)  is 
the  transmitted  wavefcm.  Although  we  cannot,  strictly  speaking, 
use  Eq.  46  in  Eq.  45  directly  because  it  is  not  a  continuous 
function-  we  cculd  replace  the  delta  function  by  an  approximating 
sequence  of  continuous  functions.  However,  the  formal  manipulations 
are  perhaps  more  clear  if  we  are  less  rigorous.  Thus  we  introduce 


1 


the  generalized  function  defined  by  Eq.  46  in  Eq.  45.  In  order 
to  carry  out  the  integration  over  the  delta  functions  in  Eq.  45 
it  is  necessary  to  revert  to  (p>x>y,z)  as  variables  of 
integration.  Since  we  do  not  in  general  have  a  one-to-one 
correspondence  between  (t,x,y,z)  and  (o,p,cp,8)  a  point  (t, 0,0,0) 
may  be  covered  by  many  points  (a»p»cp»9).  We  assume  there  are 
finitely  many.  In  addition  we  assume  that  the  origin  is  isolated 
from  caustics  of  .  Each  of  the  points  (o>p>cf>8)  covering 


(T)0,0,0)  then  has  a  neighbourhood  U  that  has  a  one-to-one 


mapping  onto  a  neighbourhood  of  (t, 0,0,0).  Then 


u(p)  =  Tj  p  p  f  p  S(t)  6(*>y,*)  Vo0-"  6(1 -p)  D(^yp  dodpdcpd0 

U  W  t  U  r  •  ^ 

u 

n  [Eq.  47; 

=  f  r  6(X)y,z)  Vo  o”“  6(1  -  p)  dt  dx  dy  dz 


n 


Now,  considering  t  to  be  a  function  of  (p,x,y,z)  we  write 

t  =  tQ  -  T(p,x,y,z)  [Eq.  48 1 


where  T  may  be  interpreted  as  the  travel  time  along  the 
geodesic  emanating  from  the  point  P.  Then 


u(P)  =  -  ^  S  jjjj’  S(t0- T)  6  (x,y,z)  V0  a  2  6(1  -  p)  Tp  dp  dx  dy  dz 


V  ' 


Eq,  49' 


where  is  the  travel  time  along  the  n-t-h  ray  from  (x0,y0,z0) 


to  the  origin,  and  nn  is  the  corresponding  value  of  o* 


From  Eqs.  18,  19  and  20  t  and  r  are  determined  as  functions  o' 
(o>p  jcp,  8)  •  We  write 


t  =  t(o,p,co,8) 


r(o  ,p  ,cp,  9) 


Eq.  50] 


Eq.  511 


Then 


21  =  -2l  _  „  T,  E 

dp  dp  dp 


[Eq.  52] 


Now  \Jj  =  t J.  T(p,x,y,z)  is,  for  fixed  p  ,  an  integral  surface 
of  the  linear  partial  differential  equation 


JL  .  2 

c5  -  (vf)2  =  constant 


fEq-  53] 


whose  characteristic  strips  are  generated  by  the  geodesic 
Eqs.  18,  19  and  20.  Comparison  with  Eq.  11  shows  that  the 
constant  on  the  right-hand  side  of  Eq.  53  must  be  (1  -  p2 )/c02 
and  viji  =v/o0  .  Hence 


c_  vT  = 


-  dr 


da 


Thus,  from  Eq.  52 


_c  21  =  c  21  +  .  2£ 

0  dp  °°  dp  dp  So 


=  c  £1  +  21.  v 

c°  dp  dp 


Now,  differentiating  Eq.  55  with  respect  to  o  we  have 


A  (  _c  21) 

do  1  c°  dp' 


c  -ill  +JlS-  .  v  +  2£  .  2i 

0  dp  So  dp  do  So  So 


Let  us  temporarily  write  (c0/c)2  =g(x,y,z).  Then  the 
geodesic  Eqs.  18 ,  19  and  20,  become 


21  = 


'°  So 


P 


=  4  vb 


So 
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fEq.  541 


[Eq.  55] 


[Eq.  56] 


[Eq.  57] 


[Eq.  53] 


From  Eq.  11  we  obtain 


P  -  v2 


1-P8 


[Eq.  59] 


v  -  p  +  p2  -  1 


Referring  again  to  Eq.  56  we  then  have 


-i-  (_  c  Bl)  = 

a a  K  0  ap; 


.  i&  +  i  1  fv2)  +  A£.  (.lnD\  = 

aP  2  aP  ^  '  +  aP  12  7  V’ 


[Eq.  60] 


Hence 


„  c  a± 

0  sp 


[Eq.  61] 


since  the  relation  clearly  is  valid  for  small  o.  Thus  Eq.  49 
becomes 


u(P)  - 


“  4r?  ^  s<*o  “  Tn^  Vo/(co  °n) 


r Eq.  62] 


where  V0  is  given  by  Eq.  38. 


The  Riemann  coordinates  involved  in  the  definition  of  J 
[Eq.  15]  are  given  in  terms  of  (g,p,cp,8)  by 


-  -cl  Cf 


=  po  cos  0  cos  cp 


=  pa  cos  9  sin  cp 


-  pa  sin  9 


[Eq.  63] 
[Eq.  64] 
[Eq.  65] 
[Eq.  66] 


Then 


J 


/t  x  y  z  > 
DW°  x1  x2  xa' 


D(txyz 
'0  p  cp  0 


)t< 


X°  X1  X2  X11 


p  cp  9 


s) 


[Eq.  67] 


■D(a^)  -  <=o/(p2o^oose) 


KV*1' <po8  003  e) 


Hence,  setting  P  —  1  , 


and  Eq.  62  becomes 


u(P) 


_1_ 

4tt 


£  S(t0  -  T  ) 
n  n 


[Eq.  68] 


[Eq.  69] 


This  is  the  geometrical  acoustics  solution  in  its  generalized 
form. 
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DISCUSSION 


The  author  confirmed  that  the  signal  distortion  can  be  obtained 
directly  if  the  source  function  is  bounded. 


B.O.  Koopman 
Arthur  D.  Little,  Inc. 
Cambridge,  Mass.,  U.S. 


The  investigation  that  I  wish  to  outline  here  originated  in  the 
military  need  of  forecasting  detection  probabilities  of  underwater 
acoustic  emitters,  in  cases  in  which  only  a  somewhat  general  idea 
of  the  sound  speed  profile  is  available,  based,  for  example,  on  a 
known  geographical  location  and  average  sound  speed  behaviour  at 
a  given  time  of  year,  Sint a  detection  of  a  distant  object  is 
the  objective,  only  very  faint  signals  enter  —  certainly  nothing 
that  could  produce  the  non-linear  effects,  shock-waves,  etc., 
that  have  been  discussed  ;n  many  of  the  papers  reported  here. 
Evidently  the  object  of  present  interest  is  not  this  or  that 
result  of  meticulously  accurate  computations  based  on  exact 
knowledge  of  the  sound  speed  c  as  a  function  of  position,  but 
more  general  facts  that  are  relatively  st able  —  i.e.,  are  not 

radically  altered  by  slight  changes  in  the  function  c.  Moreover, 
it  is  not  only  necessary  that  the  stable  evaluations  be  rather  rough 
approximations  (since  we  cannot  know  values  of  c  in  an  area  of 
future  enemy  operations  except  roughly  )  —  but  this  is  sufficient 
for  military  applications. 

These  requirements  lead  us  in  two  directions:  quantitative 
generalities:  mathematically  this  means  theorems  rather  than 

detailed  computations ;  and  a  statistical  attribute  of  the  results. 
Our  situation  suggests  a  similar  one  in  the  quantitative  study  of 
those  other  complex  physical  systems,  composed  of  the  molecules  of 


*  The  investigation  reported  here  grew  out  of  an  earlier  phase,  supported  by  ohe 
U.S.  Naval  Ordnance  System  Command.  The  author  is  uniquely  responsible  for  all 
,/iews  expressed  herein. 
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a  gas;  and  just  as  statistical  mechanics  bases  its  methods  on 
Hamiltonian  theory,  phase-space  and  its  integral  invariants, 
so  wt  shall  find,  rather  surprisingly,  that  similar  branches  of 
classical  mathematics  will  play  an  essen*  '.l  part  in  our 
investigations . 

It  is  not  too  surprising  that  the  methodology  of  classical 
Hamiltonian  dynamics  should  enter  our  problems  of  hydro-acoustic 
propagation:  mathematics  does  not  know  the  difference  between 
Fermat's  principle  of  least  time  and  Maupertuis*  principle  of 
least  action,  which  led  immediately  to  the  Hamiltonian  theory. 

But  let  us  stare  from  the  beginning. 

Our  starting  point  is  D'Alembert 1  s wave  equation  in  the  velocity  potential 

-  *tt/c2  =  o 

(with  possible  slight  modification  in  derivatives  of  lower  order) 
together  with  the  energy  density  expression 

E  “  -y-  [  |Vi)i |3  +  |ft  |7c2] 
and  the  energy  flux  vector 

*•*  =  -P|?t  **|  • 

Note  that  we  are  using  the  absolute  values,  to  allow  for  complex 
wave  functions  \|t.  We  recall  that  for  every  c  which  is  a  function 
of  geometriccl  position  and  is  independent  of  time  t,  the  wave 
equation  has  a  consequence  that  an  equation  cf  continuity  is 
satisfied;  i.e.,  that 

dE/at  +  7'?  =  0  . 

Geometrical  acoustics  is  a  valid  approximation  at  sufficiently 
high  frequencies:  U)»c/depth,  where  u>  =  2rt  X  frequency .  The 
element  which  transports  hydro-acoustic  energy  is  the  travelling 
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wave,  that  is,  a  solution  of  the  wave  equation  of  the  form 


'll  =  uf[w(t  -  S)  ] , 


S  =  S (x,y , z )  . 


Here  the  function  f  =  f(x)  must  be  defined  for  all  x.  For 
monochromatic  steady  state  propagation,  we  take  f(x)=elx,  while 
for  an  infinite  pulse,  f(x)  —  6(x),  the  Dirac  delta  function,  etc. 
The  coefficient  u  depends  on  x,  y,  z,  t,  and  even  :uj  but  it  is 
thought  of  as  varying  "slowly"  with  these  quantities.  When  the 
above  ^  is  inserted  in  the  wave  equation  and  only  terms  in  U)S 
retained,  the  eikonal  equation  is  obtained, 


|*S|3  _  l/c2  =  0  . 

A  solution  S  =  S(x,y,z)  of  this,  when  set  equal  to  a  constant 
has  for  locus  a  wave  front;  and  the  family  of  such  loci:  S(x,y,z)  =t 
is  a  moving  surface  as  the  time  t  increases —the  wave  front  of 
the  travelling  wave.  In  space-time  it  is  a  characteristic 
hypersurface  of  the  wave  equation. 

The  eikonal  equation  has,  in  its  turn,  characteristic  manifolds, 
the  bicharacteristics  of  the  wave  equation  or  rays.  The  classical 
theory  of  all  these  relations  —  known  for  well  over  a  century  — 
gives  us  the  rule  for  writing  the  differential  equations  of  the 
latter.  We  replace  3S/3x,  etc.,  in  the  eikonal  equation  by  p^, 
etc . ,  and  set 

H  -  4(p“  +  Py  +  p\  -4)  ’ 

so  that  the  eikonal  equation  could  be  written  as  H  =  0.  Then  we 
write  the  system  of  six  differential  equations 


dT=^L  = 

N*  *  tt.  r»  l)  rt  rv 


£ 

i 


5P|^iS  i?1 r£> 


where  the  denominators  are  the  partial  derivatives  of  H  with 
respect  to  the  six  independent  variables  (x,y , z, p  , p  , p  ) ,  and 

A  J  ^ 

T  is  a  parameter. 

Clearly  the  equations  for  the  bicharacteristics  are  in  Hamilton's 
canonical  form; 


dx  _  dH  , 

-5—  =  =p  ,  etc., 

dT  dp  *x 


.  |H  =  A  JL  ,  etc 

dx  5x  2C3 


They  are  the  equations  of  the  motion  of  a  particle  of  unit  mass, 
solicited  by  a  force  derived  from  the  potential  -l/2c2,  referred 
to  a  "time-  parameter  t.  We  shall  call  it  the  pseudo-particle 
and  t  the  pseudo-time . 


If  in  these  equations  the  three  "momentum"  variables  Px,Py.,Pz 
are  eliminated,  we  find  a  system  of  three  differential  equations 
of  the  second  order.  If,  finally,  we  replace  the  independent 
variable  dT  by  the  pnysical  time  t  by  means  of  the  relation 
dT  =  c2dt,  our  equations  become  identical  wiuh  the  differential 
equations  of  the  rays,  as  obtained  by  minimizing  the  time 
J1  ds/c  in  Fermat's  principle. 

Further  relationships  now  become  clear:  We  have  p  =  dx/dx  =  dx/cds 

A 

=  cosa/c,  cos  a  being  the  first  direction  cosine,  etc.  Thus  the 
"momenta"  are  directional  quantities  along  the  rays.  The  pseudo¬ 
velocity  is  seen  to  be  ds/dT  =  l/c  =  v,  the  refractive  index. 

This  is  all  part  of  the  wave-particle  complementarity ,  which  worried 
physicists  as  far  back  as  the  Newton-Huygens  arguments  about  light. 

What  is  the  relationship  between  the  rays  and  the  hydro- acoustic 
energy?  The  answer  is  given  by  going  back  to  the  expressions  for  the 
density  oF  energy  E  and  its  flux  vector  density  and  replacing 

iji  by  its  travelling  wave  expression.  On  discarding  lower  powers 
of  U)  and  then  using  the  eikonal  equation,  we  find  simple  expressions 
for  these  quantities  in  terms  of  the  intensity  juj2  ;  the  fact  that 


T-frWj  f*  (?) 


the  energy  flux  vector  ^  is  in  the  direction  of  the  momentum 

(p„>  P„j  p,.)>  that  is{  of  the  tangent  to  the  rays,  becomes 
*  y  ^ 

evident  from  this  substitution .  Finally,  the  equation  of  continuity- 

obeyed  by  E,?  leads  to  the  following  one  in  the  present  picture: 

If  the  pseudo-time  t  is  used  instead  of  t,  and  if  p  ,  p  ,  p, 

x  y  x 

are  regarded  as  the  components  of  a  fluid  velocity  field  (based  on 
the  function  S  we  an  using),  then  the  energy  (in  these  units) 
obeys  the  classical  equation  of  continuity. 


We  may  forecast  one  of  our  results  in  the  following  terms;  If  in 

the  above  picture  of  a  spatial  flow,  the  fluid  were  incompressible, 

then  the  energy  density  would  be  a  "first  integral"  of  the  ray 

equations;  i.e.,  it  would  remain  constant  along  each  ray,  so  that, 

by  tracing  it  to  its  source  (the  emitter)  where  its  value  is 

regarded  as  known,  we  would  have  its  value  at  the  point  in  space 

of  interest  (the  detector) ;  This  would  enormously  simplify  our 

Droblem.  But  since  the  above  flow  in  xyz-  space  is  very  far  from 

[*1*1 

incompressible,  the  above  method  is  totally  inapplicable1  •  However, 

by  using,  not  a  single  travelling  wave,  but  a  statistical  ensemble 

of  such  waves,  randomly  out  of  phase,  we  can  easily  establish  an 

equation  of  continuity  in  the  5-dimensional  "space"  of  values  of 

the  six  variables  (x,y,z,  p  ,p  ,p,  )  which  satisfy  the  equation 

xyz 

H-0.  Now  the  Hamiltonian  theory  comes  to  our  aid,  showing  that 
this  flow  in  5-dimensions  is  incompressible .  This  is  the 
consequence  cf  Liouville*  s  theorem,  of  fundamental  importance  in 
classical  statistical  mechanics.  Therefore  the  energy  density  is 
constant  along  each  bi-characteristic,  or  ray  in  the  5-dimension 
representation . 

The  "model"  of  the  action  of  the  ocean  in  transmitting  acoustic 
energy  over  long  ranges  from  an  emitter  of  naval  interest  to  a 
receiver,  experiencing  all  the  random  visc-issitudes  of  the 
environment  as  well  as  of  these  two  objeecs,  is  a  statistical 
ensemble  of  travelling  waves  $n,  and  only  approaches  a 

single  one  (a  point  source,  plane  wave,  etc.)  in  the  limit. 

It  happens  that  for  the  present  purpose  it  is  easier  to  deal  with  the 
ensemble  (however  near  to  its  limit)  than  the  limit  itself. 
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There  is  a  two-fold  situation  that  may  appear  paradoxical:  First, 
the  sum  of  two  travelling  waves  is  not  in  general  a  travelling 
wave  and  does  not  have  a  wave  front  in  the  usual  sense.  This  fact, 
which  should  be  evident  from  the  analytical  expressions,  has  too 
frequently  been  overlooked  and  has  led  to  errors  in  some  standard 
text  books  in  acoustics.  Second,  in  spite  of  the  inapplicability 
of  the  wave  front  picture,  the  ray  picture  and  the  Hamiltonian 
form  of  the  equations  continue  to  be  valid.  Moreover,  the  additivity 
of  energies  transmitted  along  several  intersecting  rays  is  an 
immediate  consequence  of  their  corresponding  to  travelling  waves 
that  are  randomly  out  of  phase.  If  they  were  in  phase  their 
amplitudes  would  add  vectorially:  there  would  be  interference. 

This  is  not  observed  under  the  physical  conditions  of  our  military 
situation,  with  its  multi-path  reception,  etc.,  lending  support  to 
the  validity  of  our  model. 


From  the  ensemole  we  are  led  to  the  replacement  of  tne  E 

and  t  ,  which  were  functions  of  position  only,  to  corresponding 
quantities  which  depend  on  ray  direction  (momentum)  as  well. 

Let  (dV)  be  an  element  of  volume  in  the  x  y  z  -  space  and  let 
be  a  direction  (a  point  on  the  unit  sphere).  If  (dD)  is  an 
elementary  cone  of  directions  (elementary  area  on  the  sphere) 
containing  3,  we  shall  select  the  sub-ensemble  of  l*nl'  consisting 
of  those  waves  whose  individual  wave  fronts  at  (dV)  have  a 
direction  in  (dD) .  To  be  precise,  this  has  to  be  required  at  some 
chosen  reference  point  in  (dV);  but  the  results  depend  only 
infinitesimally  on  its  exact  choice.  Consider  the  sum  of  energies 
in  (dV)  contributed  by  all  the  members  of  the  sub-ensemble  j*||  |  : 
to  quantities  of  higher  order,  it  is  proportional  to  the  volume  dV 
and  the  rea  dD,  and  may  be  written  as 


EdVdD  =  E(x,  y,  z,  Px>  Py>  P.£)  dV  dD  , 

where,  of  course,  p^,  py,  p^,  are  the  components  of  the  vector  3/c  . 
The  coefficient  E  in  this  expression  is  the  energy  density  per  unit 
5-volume  in  the  phase-space  M5,  When  integrated  over  the  whole  unit 
sphere  it  gives  the  energy  density  of  the  radiation  field  in  ordinary 


3-space.  _  | 
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A  corresponding  definition  is  given  for  the  energy  flux  density 
in  M5.  Then  by  an  additive  process,  starting  with  the  ordinary 
equations  of  continuity  for  the  individual  elements  in  our 

ensemble,  and  making  obvious  assumptions  about  the  latter,  an 
equation  of  continuity  in  the  phase  space  M5  is  obtained. 

We  give  only  one  result  at  this  point:  Let  dS  be  an  element 
of  surface  in  ordinary  3-space  and  dD  an  element  of  directions. 
If  Q  is  the  angle  between  the  normal  to  dS  and  the  direction 
3  in  dD,  the  rate  of  flow  of  energy  across  dS  due  to  the  waves 
of  ray  directions  in  dD  is  given  by 

E  cos  9  dS  dD  } 

where  E  is  the  value  calculated  at  a  reference  point  in  dS 
and  the  direction  3  .  [Note  that  any  change  in  3  and  dD 
during  the  elementary  time  dt  considered  has  only  a  higher 
order  effect  on  the  flux.]  This  result  corresponds  to  the 
following:  Normal  to  Momentum 


[d02]/2  =  cos  0  dS  dD 

c 

We  shall  find  an  invariantive  expression  for  the  cos  0  dSdD/c  , 
which  will  lead  to  the  desired  conclusion  of  the  behaviour  of 
E  along  each  ray. 

We  return  to  our  Hamiltonian  equations.  The  six  first-order 

differential  equations  determine  one  and  only  one  integral  curve 

through  Jich  point  of  the  six  dimensional  phase  space  of  the 

variables  x,  y,  z,  p  ,  p  ,  p  .  Furthermore,  as  the  independent 

x  y  z 

variable  t  increases,  each  point  moves  continuously  along  its 
integral  curve  (the  ray).  This  gives  rise  to  a  continuous  one- 
parameter  group  of  transformations  —  a  "flow".  Liouviile's 
theorem  declares  that  this  flow  is  incompressible. 
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The  only  part  of  the  six-dimensional  phase  space  is  the  sub¬ 
space  M5  of  five  dimensions  defined  by  the  equation 

M  :  H(x,y,z,  p  ,  p  ,  p  )  =  0 

x  y  z 

Since  the  Hamiltonian  H  is  a  first  integral  it  remains 
constant  along  each  integral  curve.  Thei efore  the  locus  M5 
is  invariant  in  the  flow,  and  contains  the  whole  of  every 
integral  curve  containing  a  point  in  common  with  it.  Thus 
there  is  produced  a  flow  in  the  sub-phase-space  M&.  As  a 
corollary  of  Liouville's  theorem  it  too  is  conservative. 

The  concept  of  integral  invariant  was  developed  at  the  end  of 

I*  2 1 

the  last  century  by  H.  Poincare  in  connect:  n  with  dynamical 
systems,  and  greatly  extended  and  given  wide  applications  by 
E.  Cartan1  1  in  tne  early  decades  of  this  century.  In  its 
simplest  forms  it  is  a  familiar  notion:  the  density  p  of  a 
fluid,  when  integrated  over  any  given  volume  of  the  latter, 
remains  invariant  during  the  flow,  sir.ce  the  mass  it  represents 
is  conserved.  In  a  perfect  fluid,  the  integral,  around  a  closed 
curve  drawn  in  the  fluid,  of  the  tangential  component  of  the 
velocity,  is  invariant  during  the  flow:  this  is  Lagrange's 
theorem  of  the  constancy  of  "circulation" ,  which  is  basic  to  his 
theory,  as  well  as  to  Helmholtz's  theory  of  vortices.  In  this 
case,  for  invariance,  the  integral  has  to  be  taken  about  a  closed 
curve:  it  is  called  a  "relative"  integral  invariant.  When  this 

condition  is  unnecessary,  the  invariance  is  called  "absolute". 

An  example  of  an  absolute  integral  invariant  is  obtained  by 
applying  Stokes'  theorem  to  the  circulation,  thus  expressing  the 
circulation  as  a  surface  integral  of  the  normal  component  of  the 
curl  of  the  velocity  —  the  vortic ity  (tourbillon) . 

The  Hamiltonian  theory  gives  us  a  set  of  integral  invariants, 
starting  with  che  basic  relative  one.  Setting 

0  =  Pc  6  x  +  Py  6y  +  pz  6z 
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and  using  the  notation  dQ  for  the  curl,  and  bracketed  product 
for  "outer”  or  Grassmann  product  (actually  a  determinant  opera¬ 
tion)  we  have  the  series 


|o,  ^  do,  r  [do an]  =  P  [do=i 


IfdoBT,  JudS* 


The  existence  of  th.i  six-dimensional  one  is  simply  the  statement 
of  Liouville's  theorem,  and  the  f ! ve-dimensional  one  is  a  direct 
consequence.  The  "density"  u  in  the  latter  case  can  be 
computed  as  a  simple  single-valued  expression  in  terms  of  the 
6-gradient  of  H. 

For  our  purposes  the  integrand  of  the  four-dimensional  integral 
invariant  is  of  particular  importance  because  of  its  simple 
geometrical  interpretation.  We  have,  in  fact,  when  the 
4-dimensional  element  in  M5  is  taken  as  the  pair  (dS,  dD) 
used  before,  that 


[dn2]  = 


-g-  cos  9  dS  dD 
c 


This  is  derived  from  the  expression  for  the  quantity 


£[d(f]  =  [dp  dx  dp  dyj  +  [dp  dx  dp,r  dz]  +  [dp  dy  dp  dz] 


where  each  bracket  is  the  Jacobian  determinant  cf  the  indicated 
quantities  with  respect  to  the  four  parameters  in  the 
representation  of  the  surface  element  in  question. 


Since- as  we  have  seen,  the  flux  of  energy  across  the  above 
element  is  given  by 


And  since  as  stated  before,  this  flux  the  integrand  of  an 
integral  invariant  —  by  the  ooiuum  vat ton  of  onorgy  —  it  follows 
by  the  general  theory  (actually  by  Cartan's  theorem)  that  the 
ratio,  nameiy  G  ,  in  a  scalar  Invariant,  constant  along  each 
ray.  By  referring  it  back  to  its  point  of  contact,  with  the 
emitter,  its  vaiue  can  be  determined.  By  doing  this  for  each 
ray  through  the  receiver,  the  total  energy  received  can  be  found 
hy  integration.  Since  one  usually  assumes  the  initial  values 
of  G  at  the  receiver  constant,  the  above  process  reduces  to 
that  of  finding  the  solid  angle  subtended  at  the  receiver  by  the 
directions  of  those  rays  which  connect  it  with  the  emitter. 

The  most  limited  view  of  the  above  results  is  that  we  have 
established  the  validity  of  a  ray  tracing  process  in  the  case 
of  a  general  sound  speed  function  c  c- c  ( x,  y  , /. )  1  Actually  we 

have  done  mere,  we  have  laid  the  basis  for  a  statist. leal 
treatment  of  perturbations  of  the  system.  But  this  cannot 
profitably  be  discussed  in  its  general  terms  in  this  Conference, 
we  shall  sample  it  in  a  simple  special  case  below.  Before  leaving 
this  subject,  it  is  noted  that  when  our  equations  are  w t  (ten  in 
general  curvilinear  coordinates,  our  densities  become  multiplied 
by  the  factor  jg  ,  where  g  is  the  determinant  ol  the  matrix 
of  coefficients  in  the  general  expression  for  the  length  squared, 
ds2  .  With  coordinates  appropriate  to  cylindrical,  spreading, 

^g=r2sin0,  etc.  This  in 
combination  with  the  invariance  of  the  energy  density 
automatically  introduces  the  appropriate  spreading  factors, 

1/r,  1/r2 ,  etc. 

Finally,  we  note  that  in  our  use  of  integral  invariants,  they 
are  understood  in  Cartan'  s  sense  —  ''sliding  invariants’, 
remaining  constant  when  all  points  are  slid  an  arbitrary  amount 
along  the  integral  curves,  without  requiring  to  be  moved  to 
synchronous  points,  as  requited  by  the  Point: a  i  f.  conception. 

We  turn  now,  to  illustrate  the  ideas  graphically,  to  the  special 
case  (to  which  so  much  of  present  ray-tracing  is  confined! )  in 
which  we  have  a  "fixed  pro tile”,  with  c  depending  on  depth  alone, 
c  =  c(z) .  Then,  as  is  generally  known,  the  equations  can  be 
integrated  explicitly,  requiring,  however,  a  set  of  numerical 
integrations  of  numerically  given  functions.  For  purposes  of 
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for  spherical  spreading, 
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illustration,  we  shall  exhibit  the  phenomena  graphically,  after 
the  easy  reductions  have  been  made. 

The  first  reduction  in  this  case  replaces  the  geometrical  3-space 
by  a  vertical  plane,  che  (x,z)  plane  or  more  appropriately  to 
cylindrical  spreading,  the  (r,z)  plane.  The  two  momenta  (pr>  p J 
and  equation  H  =  0  show  that  the  phase-space  becomes  a 
3-dimensional  one,  M3  instead  of  M5  .  This  makes  it  possible 
to  draw  diagrams  of  it  on  paper. 


is  the  validity  of  Snell’s 


p^  is  constant  along  each 


The  next  simplification  when  c-'c(z) 
law  in  the  large,  which  states  thal. 
ray.  This  constant,  which  we  denote  by  w  ,  is  called  the 
Snell  constant.  The  set  of  all  rays  in  M3  having  a  given  Snell 
constant  are  on  the  locus  of  the  equation 


w2  =  1/c2 


This  makes  it  convenient  to  use,  for  specifying  points  in  M3  , 
the  three  coordinates  (r,z,p  ).  This  M„  is  shown  in  Fig.  1 
drawn  with  these  three  coordinates  as  rectangular. 

This  is  not  a  tank  of  water  —  the  ocean  •»  s  represented  as  the 
(r,z)  plane  —  but,  if  a  tank  of  anything,  a  tank  of  phase  space. 
Eut  the  flow  in  M3  is  incompressible. 

Since  the  above  equation,  yiven  the  Kixed  Sx.ell  constant  value 
w,  does  not  contain  r,  its  locus  is  a  horizontal  cylindrical 
surface  whose  elements  are  parallel  to  the  axis  of  r.  Each  ray 
with  uhis  value  of  w  winds  aroui  the  cylinder  in  helix-like 
fashion. 


Let  a  plane  be  drawn  perpendicular  to  the  r-axis.  It  cuts  all 

the  rays  in  one  and  only  one  point.  As  the  value  of  r  at  its 

intersection  increases,  i.e-,  as  the  plane  is  moved  along  the 

r  axis  and  always  perpendicular  to  it.  the  points  of  intersection 

of  a  given  ray  >,,ve  in  this  plane;  such  a  plane  is  called  the 

surface  of  section,  and  was  introduced  in  the  study  of  dynamical 

systems  by  H  koine. aid,  and  lauer  by  G.D.  Birkhoft  and  the 
r  r-> 

.•.•abhor  Hie  x  '-m. sj  onr.ab  ion  i  nduced  by  the  rays  when  r  is 
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changed,  as  described  above,  can  be  pictured  as  a  flow  in  the 
lane.  During  this  flow,  the  representative  point  of  each  ray- 
moves  about  its  curve  of  given  Snell  constant.  Finally,  the  flow 
is  Incompressible.  This  is  a  consequence  of  the  sliding  integral 
invariant  ^  dQ  ,  whit  h,  when  evaluated  on  a  region  of  the 
surface  of  section,  is  equal  to  its  area. 

Since,  as  in  the  more  general  case,  energy  flux  is  a  constant 
along  the  rays,  the  application  of  the  relations  just  outlined 
can  serve  as  a  basis  for  the  study  of  propagation,  showing 
shadow  zones,  t x .  It  may  be  noted  that  caustics  have  disappeared 
in  this  representation;  they  re-appear  only  when  we  project  sets 
of  rays  of  M3  onto  the  (r,z)  plane. 

A  case  of  •’Teat  practical  importance  is  that  in  which  the  ranges 
are  very  a  rge  —  quite  a  number  of  ray  periods.  Then  the 
practical  uncertainty  of  the  exact  lengths  of  periods  leads  to 
the  replacement  of  the  energy  flux  density  by  its  average  over 
a  period.  Either  on  this  Lu_is,  or  by  reasoning  based  on 
"ergodic  mixing",  wc  are  led  to  consider  an  energy  flux  density, 
which  depends  on  the  Snell  constant  w  only,  and  therefore 
has  the  path  curves  on  the  surface  of  section  as  its  level  lines. 
With  a  slightly  higher  degree  of  perturbing  influences,  this  flux 
becomes  essentially  constant  over  those  parts  of  the  surface  of 
section  where  long-range  propagation  is  not  intercepted  (by 
underwater  obstructions,  etc.)  —  and  zero  over  the  latter  parts. 
Then  the  acoustic  power  born  by  a  bundle  of  rays  is  proportional 
to  the  (invariant)  area  in  which  it  is  cut  by  the  surface  of 
section.  On  reducing  the  picture  back  into  geogra >h Leal  space 
the  spreading  factor  comes  out  of  the  equations  automatically. 
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NOTES  AND  REFERENCES 


X.  It  would  be  sufficient  if,  instead'  of  being  incompressible, 
the  flow  in  the  xy  2-space  were  conservative,  i.e.,  had  a 
density  function  whose  integral  over  any  piece  of  this  space 
remains  constant:  then  the  ratio  of  the  energy  flux  density 
to  the  latter  would  be  the  required  first  integral.  Such  a 
density  exists  —  and  is  put  in  evidence  by  convential  ray- 
tracing  —  as  long  as  we  stay  sufficiently  close  to  the 
emitter.  This  provides  the  justification  of  the  standard 
methods;  but  orly  under  this  proviso.  Further  away,  in  fact 
at  distances  of  particular  interest,  the  densities  become 
increasingly  multiple- valued  (indeed,  singular  at  the  branch- 
loci,  the  caustics);  therefore  the  justification  breaks  down. 

It  is  for  this  reason  that  the  present  approach  is  not 
submitted  merely  as  an  alternative  to  a  more  conventional 
one,  but  as  a  method  of  salvaging  the  latter  when  it  ceases 
to  be  applicable. 

2.  H.  Poincarfe,  "M6thodes  Nouvelles  de  la  M6canique  C61bste" . 
Vol.III  (Gauthier-Villard,  Paris,  1899)- 

3.  E.  Cartan,  "Leqons  sur  les  Invariants  Intfegraux"  (Hermann, 
Paris  1922)  . 

4.  Thus  justifying  the  conventional  ray-tracing  method  in  the 
neighourhood  of  the  emitter,  and  replacing  it  by  a  method 
that  is  valid  at  greater  distances.  See  Ref.  1  above. 

5*  H.  Poincare,  l.c.j  G. D. Birkhof f ,  " The  Restricted  Problem  of 

Three  Bodies"  (keadiconti  del  Circolo  Matematico  di  Palermo, 

23  August  1914);  B.O.  Koopman,  "On  Rejection  to  Infinity 
and  Exterior  Motion  in  the  Restricted  Problem  of  Three 
Bodies",  (Trans.  Amer.  Math  Soc  .  ,  ll>27). 

6.  The  simplified  picture  (given  in  so  many  treatments)  of  the 
"limiting  critical  ray"  composed  of  *-wo  tangent  circles,  one 
above  and  one  below  the  horiaonoal  ray,  is  derived  from  the 
two-line  " aporoximation"  to  the  acoustic  profile.  Unfortunately 


for  this  picture,  the  differential  equations  determining 
the  rays  involve  the  deri vative  of  the  profile  as 
coefficients.  Since  the  two-line  11  approximation"  has  no 
derivative  at  the  point  at  which  this  critical  ray  is 
constructed,  it  is  difficult  to  understand  the  logic,  of 
the  construction . 

Instances  of  the  surface  of  section  are  shown  in  Fig.  2 
(cne  duct)  and  F'g  3  (two  ducts).  These  show  the  lines 
of  the  2-dimensional  flow,  along  each  of  which  the  Snell 
constant  w  has  a  fixed  value.  They  enclose  the  ducts. 

In  the  two  duct  case,  they  intersect  at  the  point  of  maximum 
sound  speed,  corresponding  to  the  unstable  horizontal  ray, 
approached  asymptotically  by  its  neighbours  (with  increasing 
or  with  decreasing  range.  In  Fig.  2  we  have  heavily 
shaded  the  region  at  which  the  emitter  injects  its  energy. 
Since  the  vertical  dimens'cns  of  the  emitter  are  small  in 
comparison  with  the  dept.  ,  the  region  is  a  slender  band. 

Its  horizontal  extent  is  wide  since  this  corresponds  to  the 
directions  (or  momentum  values)  at  which  it  emits  power. 

For  a  point  source,  the  band  would  shrink  up  to  a 
horizontal  line  segment. 

Figure  4  shows  the  effect  of  an  underwater  obstruction, 
with  a  key  to  the  calculation  or-.  the  left,  which  refers 
back  to  the  surface  of  section.  Time  d^es  not  permit  us  to 
go  into  details  here;  we  merely  note  that  the  fraction 
intercepted  by  the  obstruction  is  the  horizontal  interval 
through  which  the  ray  (in  rz-space)  can  f  oved  and  still 
cut  the  obstruction,  divided  by  the  ray  period. 

Figure  5  shows  a  graphical  method  of  exhibiting  the  source- 
to-duct  coupling.  The  case  shown  is  that  of  a  single  duct 
under  an  inversion  layer,  and  assumes  perfect  specular  reflec¬ 
tion  (plus  phase  randomization)  at  the  water's  surface.  One 
may  think  of  the  whole  diagram  as  reflection  in  the  latter 
surface  (method  of  images,  or  "Lloyd's  miiror'1),  whereupon 
it  resembles  the  case  of  three  ducts  separated  by  two  uns  °'le 
horizontal  rays.  Again  the  thin  horizontal  band  shows  t  ■<* 


energy  injected  by  the  emitter ,  its  heavily  shaded  part 
being  that  portion  that  survives  bottom  absorption  and 
can  undergo  long  range  propagation.  Their  ratio  might  be 
called  the  " source-to-duct  coupling  factor".  At  the 
receiver  is  drawn  a  horizontal  band  representing  what  it 
can  receive  (or,  by  reversal  of  path,  emit) .  The  dotted 
part  of  the  receiver  band  shows  the  part  of  the  emitted 
energy  that  reaches  the  receiver.  Clearly  the  fraction 
of  emitted  power  reaching  the  receiver  is  not  even 
approximately  equal  to  the  product  of  an  emitter  coupling 
factor  (determined  by  its  depth)  times  a  receiver  coupling 
factor  (determined  by  its  depth) :  given  the  acoustic  profile, 
the  fraction  in  question  is  a  function  of  two  variables  (the 
two  depl  hs) ;  but  not  a  product  of  two  functions  of  one 
variable  each. 


DISCUSSION 

The  author  stated  that  the  theory  %vould  have  to  be  reworked 
for  the  5-dimensional  case  when  applied  to  a  rangt -dependent 
sound-speed  profile.  Also,  since  a  flat  bottom  was  assumed 
in  the  theory,  it  would  have  to  be  modified  if  that  were  not 
the  case. 
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RAYS  AND  STATISTICAL  DIFFRACTION  THEORY 


by 

R.H.  Clarke 

SACLANT  ASW  Research  Centre 
La  Spezia,  Italy 


ABSTRACT 

A  method  of  extending  ray  tracing  is  proposed,  such  that  medium- 
scale  irregularitios  of  a  statistical  nature  arc  taken  into  account 
along  with  large-scale  irregularities  in  refractive  index  of  the 
medium. 

INTRODUCTION 

Ray  tracing  is  a  relatively  simple  and  very  practical  method  of 
obtaining  solutions  to  the  fundamentally  very  difficult  problem  of 
wave  propagation  in  an  inhomogeneous  medium.  These  solutions  are 
admittedly  approximate,  but  give  perfectly  satisfactory  "engineering" 
answers,  provided  the  scale  size  of  the  inhomogeneities  is  large 
compared  to  the  propagating  wavelength,  and  provided  one  does  not 
require  answers  in  regions  of  focusing  or  shadowing. 

I  want  to  address  myself  to  the  problem  of  dealing  with  irregularitic 
both  in  the  medium  and  on  the  boundaries,  whose  scale  size  is  not 
large  compared  to  the  propagating  wavelength.  In  particular,  I  am 
interested  in  scale  sizes  of  the  same  order  or  somewhat  larger  than 
the  propagating  wavelength.  I  shall  exclude  from  consideration  those 
irregularities  whose  scale  size  is  smaller  than  the  propagating 
wavelength,  since  in  a  sense  these  will  be  invisible  to  the 
propagating  *ave  and  are  likely  to  have  only  a  collective  effect. 
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|  (For  example,  bubbles  will  have  the  effect  cf  changing  the  acoustic 

A 

l  properties  of  the  water  when  encountered  by  metre  wavelengths,  and 

|  so  the  problem  reverts  to  that  of  the  effect  of  conglomerates  of 

Yc 
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|  such  bubbles  of  scale  sizes  of  a  wavelength  or  larger.) 
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One  can  state  the  problem  in  a  somewhat  more  restricted  fashion  by 
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e  asking  how  one  can  retain  the  advantages  of  using  ray  tracing 

1 

|  to  deal  with  large-scale  irregularities  when  medium-scale 

i 
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|  irregularities  are  also  present.  I  eventually  want  to  discuss 

| 
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1  the  effect  of  medium-scale  irregularities  on  the  sea  surface, 
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|  in  the  depth  of  the  thermocline  and  in  the  medium  everywhere.  But 
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I  it  is  useful  bo  sta^t  with  the  rough  surface. 
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|  EXAMPLE  OF  SURFACE  ROUGHNESS 

i: 

|  After  a  more-or-less  ^ortuous  path,  some  of  the  rays  leaving  the 

|  source  will  '  i  fke  the  surface  [Fig.  1].  A  straightforward 

1  extension  of  the  ray  tracing  concept  would  be  to  say  that  each  ray 

|  is  reflected  in  the  local  specular  direction.  But  there  are  two 

i  objections  to  this  course:  (1)  this  idea  is  only  valid  for 

1 

| 

1 

It  large-scale  irregularities  and  (2)  a  hideously  large  nui:.ber  of 

ray  tracings  would  have  to  be  made  in  order  to  obtain  a  satisfactory 

R.  statisticai.  ensemble. 
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|  An  alternative  approach,  and  the  one  I  shall  advocate,  is  to  use  ray 

i 

F  tracing  to  just  below  the  surface;  then  to  employ  statistical 

F  # 

diffraction  theory  to  account  for  the  effect  cf  the  statistically 

% 

f  rough  surface;  and  then  to  employ  ray  tracing  again  to  describe 

i 

I 

%  the  suboequenb  progress  of  the  field.  This  gets  over  the  two 

i 

|  disadvantages  I  mentioned  in  connection  with  the  first  approach: 
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1  (1)  since  diffraction  is  taken  into  account  the  validity  of  the 

i 

*  method  is  p*»t  restricted  tc  large-scale  irregularities,  and  (2)  the 

i 

|  ensemble  averaging  is  performed  at  the  surface  and  so  a  single  ray 

j  tracing  suffices  to  describe  the  subsequent  reflected  field. 
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I  shall  give  more  details  of  this  in  a  moment.  But  first,  I  must 
describe  more  clearly  how  one  can  make  the  transition  from  rays 
to  fields  and  back  again. 

RAYS  AND  PLANE  WAVES 

.  .  Ill.l  '-IIMI..I  —  I  . .  ■  I . . — 

Dating  from  Rayleigh’s  treatment  of  the  problem  of  reflection  from 
a  corrugated  surface  [Ref,  1],  the  expansion  of  acoustic  fields 
in  terms  of  plane  waves  travelling  in  different  directions  has 
become  increasingly  popular.  Intuitively,  one  expects  there  to 
be  an  equivalence  between  such  plane  waves  and  the  purely 
geometrical  concept  of  rays; that  equivalence  will  now  be 
demonstrated. 

Uniform  Medjum 

The  pressure  field  p(x,y,z)  in  the  half  space  z>0  can  be 
represented  [Ref.  2]  by  the  plane-wave  spectrum  F(a,  g)  ,  such 
that 

CO 

PP  . 

p (x, y ,  z)  ”  jjF(a,p)  exp  l-jk0(ax+  gy  +  yz)  1  da  d8  [Eq.  1] 

—  Of' 

where  (a,g,y)  are  the  cosines  of  the  angles  formed  between  the 
direction  of  a  single  plane-wave  component  and  the  three  rectangular 
coordinate  axes  (x,y,z),  and  k0  is  the  phase  constant  (wave- 
number)  of  the  medium. 

Assuming  the  acoustic  source  to  be  at  the  origin,  at  a  large 
distance  r  from  the  source,  such  that  k0r>>l,  it  can  be  shown 
by  applving  stationary  phase  methods  that  the  pressure  field  is 
asymptotica.1  ly 

P  3^  F(«,P)  *"jkr  • 
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[Eq.  2] 


Thus  the  angular  plane-wave  spectrum  F(gj,P)  is  proportional  to 
the  directivity  pattern  of  the  source. 


The  associated  intensity  is 
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where  Z0  is  the  characteristic  impedance  of  the  medium.  It 
will  be  useful  later  to  refer  intensities  to  the  "intensity  at 
unit  distance",  which  is 
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[Eq.  3] 


Layered-Inhomogeneous  Medium 

When  the  acoustic  properties  of  the  medium  change  with  z,  then 
the  angular  spectrum  F(o:,p),  which  describes  the  field  at  the  source 
level  z  —  0,  no  longer  describes  the  field  for  any  other  z. 

However,  by  writing  the  angular  spectrum  as  a  function  of  z, 
namely  F(a,p,z),  it  can  be  supposed  that  all  plane  waves  emanating 
from  the  source  follow  the  ray  paths  prescribed  by  geometrical 
acoustics,  provided  the  irregularities  in  the  medium  are  of  scale- 
size  very  large  compared  to  the  acoustic  wavelength.  Hence  [Ref.  3], 

CO 
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P <.x,y,  z)  -  jj  F (a ,  P,  z)  exp  l- jk0(ax  +  Py  +  qdz)}dadg  [£q.  4] 

—  co  0 

where  q  =  ncos9,  the  refractive  index  n  =  n(z)  =  r(0)/c(z)  is  the 
ratio  of  the  sound  speeds,  and  0  is  the  local  angle  to  the  vertical 
of  the  ray  path  LFig.  2],  From  Snell's  law  it  is  obvious  that 


2  _  2  .2 

q  -  n  -  a  -  P  . 


[Eq.  5] 


The  W.K.B.  solution  [Ref.  4]  for  a  plane  wave  travelling  in  a 
layered-inhomogeneous  medium  yields 


F(a,  0,  z) 


Pc  cos  0O  \  lr 


P0c0  cos  0 


F(a,p) 


[Eq.  6] 


where  p=p(z)  is  the  density,  and  the  subscript  0  refers  to 
the  source  level,  (Note  that  Z  —  Pc  is  the  characteristic 
impedance  at  the  level  z) . 


Integrating  Eq .  4  by  stationary  phase  methods,  asymptotically 
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exp  {-jk0(a  0x+  0oy  -k  q  dz)  | 


[Eq.  7] 


where  (a  0, $0)  are  the  direction  cosines  at  the  source  which 
satisfy  the  stationary-phase  conditions.  These  conditions  are  that 

x  =  -  J  [Eq.  8a] 

o 

and  that 

y  =  -  P  dz  [£q*  8b] 

which  are  just  the  equations  of  the  ray  path.  The  quantity  A  is 
a  determinant  in  the  general  case  [Ref,  3]>  but  in  the  x-z  plane 
reduces  to 


A  = 
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i_adz)  (  ,  3-adz)  . 

Sa  J0 


[Eq.  9] 


The  intensity  corresponding  to  the  pressure  of  Eq.  7  is 
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[Eq.  10] 


where,  as  can  be  shown  by  applying  Eq.  9  bo  Eq.  5, 


n  z 


sin2  8ft  «■*  cos'5  9 


sin  8 


dz  . 


[Eq.  11] 


Then,  the  final  intensity  formula  is 

i  =  3o  - - 
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sm  8; 


x  cos  8  cos  80  P  'S—n ^ —  dz 


[Eq.  12] 
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and  is  equivalent  to  the  formula  developed  by  Krcl  [see  Session  2 
of  these  Proceedings]  using  purely  geometrical  arguments. 

Thus  an  equivalence  between  a  ray  description  an  a  plane-wave 
description  of  acoustic  fields  in  a  layered  inhomc  -eneous  medium 
has  been  established,  which  takes  care  of  the  effect  of  large-scale 
irregularities . 

The  next  step  in  the  argument  is  to  :onsiuer  how  each  of  these 
plane-wave  components,  which  r,o  to  make  up  the  total  field,  are 
affected  by  medium-scale  statistical  irregularities  encountered 
either  in  the  medium  or  at  its  boundaries. 

STATISTICAL  DIFFRACTION  THEORY 

Consider  the  simplest  case,  shown  in  Fig.  3>  of  a  plane  wave  incident 
normally  on  a  "random  phase  screen".  Such  a  screen  alters  the 
phase  of  a  wave  propagating  through  it  in  a  random  manner,  but 
leaves  its  amplitude  unchanged.  (The  physical  mechanisms  in  the 
ocean  which  produce  such  random  phase  screens  will  be  discussed 
in  the  next  section.) 

If  the  random  phase  indjced  by  the  screen  is  a  zero-mean,  gaussian 

random  process  of  variance  o',,  chen  the  transmitted  field  will 
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consist  of  a  cohf^snt  part  and  an  incoherent  part.  (Coherence  is 
used  here  in  the  sense  of  the  phase  having  a  deterministic  relation 
to  the  incident  phase.)  It  can  be  shown  [Ref.  5]  that  the  coherent 
part  of  the  transmitted  field  is  a  plane  wave,  in  all  respects 
the  same  as  the  incident  field,  except  that  its  amplitude  is  reduced 
by  exp{-7r  This  can  be  expressed  by  saying  there  is  a 

"coherence  loss"  of  intensity  of 

exp  1-0  or  4*34°  |  dB  . 

But  this  is  not  a  real,  absorptive  loss,  and  the  remaining 
transmitted  energy  is  incoherently  scattered  in  a  pattern  which  is 
determined  by  the  second-order  ^i.e.,  lateral  correlation) 
statistics  of  the  phase  across  the  screen. 

In  terms  of  rays:  the  incident  ray  suffers  a  "loss",  but  apart 
from  that  continues  as  though  the  screen  were  not  there.  The  lost 
energy  is  converted  at  the  screen  into  new,  incoherent  sources  of 
energy  whose  angular  pattern  can  be  determined.  Ray  tracing  can 
be  applied  to  follow  the  subsequent  behaviour  of  this  new  -ource  of 
acoustic  energy. 

RANDOM  MECHANISMS  IN  THE  SEA 
Rough  Sea  Surface 

For  a  plane  wave  incident  obliquely  on  a  randomly  rough  sea  surface, 
[Fig.  4],  the  simplest  (and  most  common)  approach  is  to  ignore 
amplitude  effects  and  to  consider  only  the  random  phase  induced  in 
the  incident  wave  arising  from  the  local  excess  path  travelled  by 
the  wave  to  and  from  the  surface,  compared  with  reflection  from  the 
mean  surface.  Thus  the  surface  is  replaced  by  a  random  phase  screen 
If  the  surface  profile  is  a  zero-mean,  gaussian  random  process  of 
variance  0^  the  random  phase  variance  is 

a  2  2 

.  0$  =  (2k  cos  8)  c£  . 
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Thus  tao  incident  ray  is  specularly  reflected  with  a  coherence 
less  of  4*34  a | ,  and  the  remaining  energy  is  scattered  incoherently 
with  an  intensity  pattern  determined  by  the  spatial  correlation 
function  of  the  surface  roughness. 

Internal  Waves 

Figure  5  shows  an  idealized  model  of  an  abrupt  thermocline 
boundary  separating  two  regions  of  the  ocean  in  which  the  sound 
velocities,  and  hence  the  phase  con  bants,  ki  and  kz,  are  slightly 
different.  If  the  boundary  profile  is  a  zero-mean,  gaussian  random 
process  of  variance  o3,  then  the  same  sort  of  arguments  used 
for  the  rough  surface  establish  the  first-order  effect  on  an 
obliquely  incident  plane  wave  of  such  boundary  as  a  random 
phase  screen  with  phase  variance 

=  (ks-ki)2  sec3  0  o3  . 

(A  similar  expression  has  been  used  to  examine  the  effect  of 
irregularities  in  dielectric  holograms  [Ref.  6].) 

Volume  Irregularities 

If  a  plane  wave  travels  a  distance  t  [see  Fig.  6 ]  through  a  slab 
of  tenuous  irregularities  in  refractive  index,  then  it  is 
physically  plausible  to  suppose  that;  the  emerging  field  is  randomly 
modulated  in  phase  but  unaltered  in  amplitude.  (For  a  more 
rigorous  validation  of  this  approach,  see  the  r6sum6  of  the  work 
of  Fejer  and  Bramley  in  Ref.  5*)  Hence  the  slab  of  irregularities 
behaves  as  a  random  phase  screen.  If  l  is  many  times  a  typical 
scale,  size,  C0,of  the  irregularities,  then  a  crude  application  of 
the  Central  Limit  Theorem  establishes  that  the  emerging  phase 
is  approximately  gaussian,  with  variance 
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where  o  is  the  variance  of  the  refractive  index  fluctuations, 
n 

Hence  as  a  ray  traverses  such  irregularities  it  will  suffer  a 
loss  of 

4.34k2  C0<^  dB/unit  length 

of  its  energy  to  incoherent  scatter,  the  angular  spread  of  which 
will  be  determined  —  as  in  the  other  examples  —  by  the  lateral 
scale  size  of  the  irregularities. 
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APPROXIMATE  METHODS  FOR  RAY  TRACING 

by 

M.J.  Haintith 

Admiralty  Underwater  Weapons  Establishment 
Portland,  Dorset,  U.K. 

INTRODUCTION 

The  conventional  approach  to  ray  tracing  is  to  follow  one  ray,  usually 
specified  by  its  initial  direction,  by  standard  techniques  of  numerical 
integration  along  the  ray  path,  building  up  a  set  of  values  of 
horizontal  displacement,  direction,  and  travel  time,  as  a  function  of 
vertical  displacement  (horizontal  stratification  of  the  sound  speed 
profile  in  the  medium  is  assumed) .  This  process  demands  computer  time 
and  storage.  By  this  means  a  family  of  ray  plots  may  be  built  up 
[e.g.  Fig.  1] 

In  many  applications,  however,  this  forward  computation  is  inconvenient, 
in  that  an  inverse  problem  requires  solution.  Examples  are: 

a.  Given  the  terminal  points  of  the  ray,  what  is  the 
grazing  angle  at  one  point  on  it  (this  frequently  occurs  in 
experimental  determinations  of  bottom  reflectivity). 

b.  Given  the  measured  travel  time  from  surface  to  bottom, 
what  is  the  true  slant  range  (e.g.  the  use  of  bottom  transponders 
in  some  navigational  systems) . 

c.  Given  the  known  slant  range,  what  is  the  true  travel 
time  (the  converse  of  (a),  also  often  encountered  in  bottom  studies). 

A  common  feature  of  these  problems  is  that  they  do  not  involve  rays 
having  turning  points  (i.e.,  the  slope  is  always  of  the  same  sign), 
and  basically  the  theory  to  be  described  is  restricted  to  this  special, 
but  important,  situation.  It  is  possible  to  extend  the  treatment  to 


rays  having  a  turning  point,  but  the  advantages  of  the  approach  are 
not  as  marked;  and  thix  extension  will  not  be  discussed  here. 

THEORY 

Figure  2  shows  the  geometry  of  the  situation.  Horizontal  range  is 
denoted  by  x,  vertical  depth  by  z,  and  the  grazing  angle  at  any 
point  on  the  ray  by  Q  .  The  terminal  points  of  the  ray  will  be 
denoted  by  (0,0)  and  (X,Z).  The  sound  speed  profile  (horizontally 
stratified)  is  assumed  known,  the  sound  speed  at  depth  z  being  c. 
The  slope  of  the  slant  range  line  80 ,  and  the  slant  range  is  D. 

The  basic  ray-tracing  equations  .'re  that  along  the  ray: 


tan  0  =  dz/dx 

cos  0/c  =  p  =  constant  (Snell's  Law) 
dz/dt  =  c  sin  0,  where  t  =  travel  time 


On  integrating  along  the  whole  path,  we  obtain 


[Eq-  i] 

[Eq.  2] 
[Eq.  3] 


X  = 


t  = 


cot  0  dz 


2 

J1  (cosec  0/c)  dz  . 


[Eq.  4] 


[Eq.  5] 


If  we  use  Eq.  2  to  express  Eq.  4  in  terms  of  c  and  p,  we  note 
that,  if  given  X  and  Z,  Eq.  4  becomes  essentially  an  integral 
equation  for  determining  p.  The  method  to  be  described  is  based  on 
noting  that  integration  is  an  averaging  process,  and  that  this 
suggests  that  we  are  i..  effect  computing  some  rather  complex 
weighted  average  of  c. 

As  far  as  Eq.  5  is  concerned,  we  may  note  that  another  average 
value  of  c  is  defined  by  travel  along  the  straight  line  path; 
since  by  Fermat1 s  principle  the  true  travel  time  represents  a 
stationary  value,  the  change  due  to  moving  to  the  displaced  straight 
line  path  can  differ  only  be  second  order  quantities. 


We  therefore  express  c  in  terms  of  its  deviation  from  the  mean 
value  over  depth  c  ,  where  c  =  ~  J*  Zc  dz  ,  and  write 
c  —  c[l + e (z)],  where  for  real  profiles 

c(z)  «  1  [Eq. 

It  is  also  convenient  to  replace  the  Snell's  law  constant  p  by 
an  angle  9  ,  defined  by  the  equation 

cos  0/5  =  p  .  [Eq. 

Since  c  is  a  value  which  actually  occurs  on  the  sound  speed 
profile,  0  is  a  real  angle  for  any  real  ray. 

The  mean  of  e  over  z  is  clearly  zero,  and  we  may  define  higher 
moments  by  such  equations  as 

e2  =  ^  e3  dz  .  [Eq. 

On  making  the  appropriate  substitutions  in  Eqs.  4  and  5,  we  obtain 

x/z  =  cote0  =  Jz(i+e)  [l-cot2  5  (2  c  +  e3)]“^dz  . 

°  [Eq. 

t  =  (cosec  5/c)  PZ  (l+e)-1  [  1  —  cot2 0  (2e  +  e3)]""2  dz  . 

Jo  [Eq. 

Equations  9  and  10  may  now  be  expanded  as  binomial  series  in  e  > 
the  results  being 

cot  0O  =  cot§[l  +  3A-  P-  cosec3  (!)  cot^  +  oCe2)]  [Eq. 

ct/Z  =  cosec  0  [l  +  P  (1  -  £  cot2  S'  +  3/s  cot*  5  )  +  oCe3)] 

[Eq. 


the  first-order  terms  vanishing  identically. 


If  we  retain  only  terms  to  the  second  order* Eqs,  11  and  12  are  very 
easy  to  invert  or  otherwise  manipulate >  with  the  following  results 
(noting  that  D  =  Z  cosec  0O  )  : 


cos  (3  =  cos  e0  (1  -  a/a  3s  cot2  80  ) 

D  =  ct  [1  +  £  Is*  {  (c2  t2  / z2  )-  3 }]  , 

and  the  inverse  cf  Eq.  14 

St  -  D[l--;'eij(D2/Z2)  _3|]  .  [Eq.  15] 

These  equations  clearly  give  <.  very  simple  answer  tn  the  problems 
cited  in  th"  introduction.  They  are  ■’iasy  to  compute,  and  require 
computer  st'  'age  'or  jnly  two  environmental  parameters,  c  and  , 
both  of  whir  re  ea.,i2,'  computed  one-  for  all  for  any  given  sound 
speed  profit- 


[Eq.  13] 
[Eq.  14] 


ACCURACY 

Equations  13  to  15  are  approximations  in  which  terms  in  e’3’  and 
higher  moments  h&i  a  been  ignored,  and  it  is  obviously  necessary  to 
determine  the  errors  /rtroduced  (and  indeed  even  to  decide  if  the 
series  is  convergent) . 

This  problem  has  been  solved  as  follows.  Considering  all  possible 
sound  speed  pre files  for  which  c  and  IF"  are  specified,  and  for 
given  values  of  X  and  Z,  for  which  of  these  profiles  will  the- 
values  of  p  or  of  t  V:iven  by  Eqs.  2,  4  and  5  have  extremal 
values?  This  is  a  variational  problem,  which  can  be  handled  by  the 
technique  of  using  Lagrangian  multipliers  for  the  equations  of 
condition.  The  result,  for  both  p  and  t  ,  is  that  extremal  values 
will  be  attained  when  c(z)  is  a  function  of  {z)  which  can  take 
only  two  discrete  values,  i.e.,  when  the  sound  speed  profile  is  that 
of  a  two-layered  environment. 


This,  however,  is  not  sufficient  to  determine  a  true  maximum,  since 
the  two-layer  profile  is  specified  by  only  two  conditions,  but  has 
three  degrees  of  freedom.  It  is  necessary  to  find  a  third  constraint, 
and  an  obvious  one  is  given  by  the  observation  that  any  real  profile 
has  bounded  values  of  c,  that  is,  that  it  has  a  maximum  and  a 
minimum  value  for  sound  speed.  It  can  now  be  shown  that  the  extremal 
values  in  this  situation  will  be  given  when  one  of  the  layers  is 
allocated  either  the  greatest  or  the  least  value  of  c,  denoted  by 


c  and  c  . 

max  mm 


[Fig.  3]. 
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From  these  extreme  profiles  it  can  be  shown  that  the  series  expansion 
is  absolutely  convergent,  and  that  the  values  of  cos  S’  and  1)  lie 
within  the  bounds  given  by  the  following  expression: 


cos  8  —  cos80  [l-^a'P  cot2  80  (1+Yj)] 

D  =  ct  [1+  |  ?  |(c2  t3  /z3  )  -  3  \  (1  +  ys  )  ] 


■»)  cot2  80  <  Ya  <  ( 


)  cot2  9, 


c  =  c  (1  +  a) 

max  v  ’ 


°min  ”  SU-b) 


It  is  clear  from  these  expression  that  the  error  is  greatest  at  the 
maximum  range,  and  falls  off  roughly  as  the  fourth  power  of  range. 

If  we  make  some  simplifying  assumptions  (basically  that  gradients  are 
never  very  large,  so  that  the  sound  speed  profile  moves  relatively 
smoothly  between  its  extremes),  it  can  be  shown  that,  to  a  reasonable 
degree  of  accuracy,  the  error  at  maximum  range  is  approximately  equal 
to  the  correction  introduced  by  adding  the  term  in  e5"  for  the 
Snell's  law  constant  cos  5,  and  is  half  the  corresponding  correction 
for  the  slant  range  determination. 


ILLUSTRATIONS 

To  demonstrate  the  sort  of  accuracy  that  the  approximations  can  give, 
a  comparison  has  been  made  between  the  results  of  an  exact  computa¬ 
tion,  using  a  digital  computer,  and  the  approximations  given  above, 
for  two  profiles  (chosen  basically  to  ease  the  digital  computer' s 
task* ) 

Figure  4  outlines  the  profiles  used.  That  marked  'typical'  has 
parameters  not  unlike  those  found  in  the  real  ocean;  the  'extreme' 
profile  was  designed  to  have  wide  limits  (a  10$  variation  in  sound 
speed)  and  incorporates  a  marked  inversion  layer. 
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Figures  5  and  o  illustrate  the  results  for  travel  time.  The  computed 
slant  range  errors  are  shown,  for 


a.  the  very  simple  formula  D  ~  ct  and 

b.  the  second-order  expression  [Eq-  ±5].  The  computer 
upper  and  lower  bounds  are  also  shown.  It  will  be  seen  that,  with 
one  second  order  correction  terms  the  error  is  at  most  9'm  in  28.9  km 
for  the  'typical  profile,  and  is  only  23m  xn  19.3km  for  the  'extreme' 
profile. 

A  similar  analysis  was  carried  out  to  compute  the  error  in  initial 
grazing  angle  as  deduced  from  th"!  Snell's  law  parameter  cos  5.  A 
summary  of  th^  results  is  given  in  the  following  table. 


TABLE  1 

ERRORS  IN  GRAZING  ANC' E 


Profile 

Range 

True  Grazing 

Error  in 

Maximum  Error 

(km) 

Angle 

Angle 

Bounds 

(deg) 

(deg) 

(deg) 

15 

15 

0.01 

±0.05 

20 

8 

0.1 

±0.17 

' Typical : 

24.5 

4 

0.27 

±0.4 

29-3 

0 

0.74 

±0.87 

10 

17.5 

0.01 

±0.15 

15 

7 

0.1 

±0.8 

' Extreme' 

16.5 

4 

0.17 

±1.2 

18.7 

0 

0.29 

±1.9 

It  will  be  noted  that  the  errors  in  grazing  angle  are  larger  for 
the  'typical'  profile  than  for  the  ' extreme' j  this  is  because  the 
horizontal  ranges  with  the  'typical*  profile  are  much  greater  than 
for  the  ' extreme' ,  and  the  strong  range  dependence  outweighs  the 
smaller  variation  in  the  sound  speed.  Even  so,  the  errors  are 
remarkably  small  over  most  of  the  range,  and  the  accuracy  everywhere 
is  probably  greater  than  is  warranted  by  the  reliability  of  the 
input  data. 
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DISCUSSION  OF  RESULTS 

It  is  apparent  that  for  most  purposes  the  errors  introduced  by  the 
use  of  this  approach  are  far  smaller  than  the  quality  of  the  input 
data  would  justify,  and  the  saving  in  computer  size  required  is 
considerable.  Furthermore  it  is  clear  that,  because  the  sound  speed 
enters  only  in  the  form  of  statistical  averages  it  is  easy  to  assess 
the  precision  to  which  individual  measurements  should  be  made. 

Again,  from  this  analysis,  it  is  evident  that  the  effect  of 
irregularities  in  the  profile  will  not  in  general  be  of  great 
importance;  this  is  a  deduc  ion  that  would  be  difficult  to  make 
by  conventional  ray -tracing  methods. 

At  first  sight  the  high  accuracy  of  this  very  simple  approximation 
seems  surprising.  The  following  argument  gives  an  explanation  for 
this  result.  In  the  integration  over  z  for  X  and  t  (Eqs.  4  and  5), 
the  order  in  which  successive  increments  are  added  is  immaterial,  and 
the  profile  can  be  redrawn  so  that  c  is  a  monotonically  increasing 
function  of  z  (this  is  the  same  as  forming  a  Lebesgue  integral). 

The  approximation  then  consists,  in  effect,  of  replacing  this 
' regularised'  profilt  by  the  constant  gradient  profile  of  best  fit 
by  least  squares.  The  shape  of  the  ray-path  will  be  quite  different, 
but  the  Snell's  law  constant  and  the  travel  time  will  be  nearly 
unchanged.  This  argument  also  shows  immediately  why  the  two-layer 
profile  gives  the  extreme  bounds,  since  this  is  the  one  which  is 
least  well  fitted  by  a  single  straight  line. 


m 


The  method  is  clearly  capable  of  extension.  For  example,  if  in  say, 
a  side-scan  sonar  the  launching  grazing  angle  and  the  travel  time  are 
simultaneously  recorded  it  is  possible  to  estimate  the  height  of  a 
projection  above  the  sea-bed  (since  in  effect  9  and  t  are  given) 
by  suitable  inversion  of  the  equations.  Such  applications  will  be 
reported  separately. 
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BOUNDING  PROFILES 


FIG.  4 
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FIG.  5 


ERROR  IN  SLANT  RANGE  DETERMINATION 


FIG.  6 


ERROR  IN  SLANT  RANGE  DETERMINATION 
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CONSIDERATIONS  ON  NUMERICAL  AND  EXPERIMENTAL  PROPAGATION  MODELS 
FOR  TWO-DIMENSIONAL  VARIATION  OF  MEDIUM  PROPERTIES 


by 

W.  Sluyterman  van  Langeweyde 

Forschungsanstalt  der  Bundeswehr  fur  Wasserschall  und  Geophysik 

Kiel,  Germany 


ABSTRACT 


One  of  the  problems  the  FWG  is  dealing  with  is  the  prediction  of  sound 
propagation  loss  in  shallow  water.  A  part  of  this  problem  is  being 
treated  by  bringing  measurements  and  numerical  computation  into 
relationship. 


First  a  formal  description  was  given  of  the  general  problem  and  then 
the  concept  of  a  partial  problem.  Some  sets  and  a  function  were 
defined: 


C  ;  the  set  of  oceanographic  conditions. 

E  :  the  set  of  acoustical  excitations  under  water. 

p  =  Cx  E  :  the  caroesian  product  of  C  and  E,  so  that 
p  =  (c ,  e)  sP  with  c  e  C  and  e  e  E 
called  the  set  of  channel  conditions. 

A  :  the  set  of  acoustical  fields  under  water 

f  :  P  •*  A  :  A  function  of  the  set  P  into  the  set  A, 
so  that  a  =  f  (p)  for  a  e  A,  p  e  P 


The  main  problem,  the  prediction  of  a  e  A  was  divided  into  two 
subproblems: 

(1)  Subproblem:  Prediction  of  p  e  P 

(2)  Subprohlem:  To  find  f  for  all  peP 

The  paper  was  confined  to  the  treatment  of  subproblem  two. 


The  methods  used  to  approach  the  problem,  relevant  to  propagation  loss 
in  shallow  water,  are: 


(1)  Simultaneous  measurements  of  oceanographic  and 
acoustical  data  in  special  areas  cf  the  North  Sea  and  the  Baltic. 

(2)  Acoustical  measurement  in  a  model  basin  with  definite 
physical  conditions  as  far  as  possible. 

(3)  Numerical  model  computations  on  a  digital  computer. 


These  three  methods  have  to  be  related  to  each  other.  The  corresponding 
sets  and  functions  are  written  C* ,  E' ,  P' ,  A'  and  f'  for  method  two 
and  C" ,  E" ,  P" ,  A"  and  f"  for  method  three. 

Taking  measurements  according  to  method  one,  finite  subsets  of  C,E,A 
are  found,  which  are  marked  by  a  bar. 

“  £ 

C  -  {Ci  li  el]  c  c 

E  =  lej|j  e  J|cE 

P  -  ipk.|keKlcCxE=i(ci,  ej)|iel,  jeJjcP 

£  =  iaklk  eKlcA 


> 


I, J,K  :  sets  of 
Subscripts 


The  function  thus  found,  f :  P  -*  A  is  a  restriction  of  f  ,  so  that 
£  -f  lor  pe  F  . 

The  next  aim  is  to  define  a  numerical  model  with  an  input  set  P" 
equivalent  to  F,  an  algorithm  f"  so  that  the  output  set  5"  is 
equivalent  and  comparable  to  A. 

The  numerical  model  (P" ,  f" ,  5")  will  be  called  an  isomorph  image 
onto  the  triple  (P,  ? ,  K)  with  the  isomorphism  (hi  )  ,  so  that 
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hj  :  p  ^.4*1  *»  p"  is  a  bijective  function  P  onto  P1' 

h2  :  X  ^ -U  »  X"  is  a  bijective  function  X  onto  X"  and 

[Remark:  g  .*  X  2i*U.  y  <$=$> 

Y  =  g(x)  for  all  y  Y,  x  e  X  and 
g(*i  )  *-=  g(x2  )  =£-  Xj  =  x?  ] 

f"[hi(p)]  for  all  p  e  P 

III 

3?" 

u 

X!* 

The  final  goal  is  to  extend  this  isomorphism  onto  whole  sets  P 
and  A  to  come  from  P  to  A,  without  knowing  f,  via  hx  ,  f" 
and  hg  1  (reverse  function  h2 ) ;  however,  in  order  to  keep 
the:  problem  tractable,  the  sets  have  to  be  decomposed  into 
further  cartesian  factors,  and  these  factor  sets  have  to  be 
reduced  by  classifications  and  by  means  of  statistical  descriptions* 

By  a  simple  concrete  example  of  6  series  of  model  basin  measure¬ 
ments  the  decomposition  of  the  set  P'  was  demonstrated  to  the 
point  that  the  internal  structure  of  the  first  function  Rx'  of  the 
aspired  isomorphism  could  be  seen  and  the  channel  conditions 
described. 

In  a  few  words,  the  ray  tracing  program  used  and  the  representa¬ 
tion  of  the  function  ?"  was  described.  The  program  traces  rays 
in  a  plane,  which  is  given  the  range-depth  coordinates.  A 
variability  of  sound  velocity  in  both  directions  is  taken  into 
account  in  a  manner  that  this  plane  is  divided  into  triangles 
with  constant  sound  velocity  gradients  with  vertical  and  horizontal 
components.  This  causes  the  rays  to  be  triangle-wise  parts  of 
circle  arcs. 


h8[f(p)]  = 
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The  equation 


V(x,2)  =  U  (x-xc)  +|2  (E-rc) 

with  x  ,  z  as  centre  coordinates  of  the  ray  curvatures,  is 
c  c 

generally  valid  for  two  dimensional  variability  sound  velocity. 

For  v  =  ax+gz+y,  x  and  z  are  constant,  which  leads 

c  c 

to  circles. 

For  the  surface  and  bottom  reflections  exits  for  subroutines  are 
provided  which  compute  for  instance  random  variable  reflection 
angles  at  surface  or  reflection  factors  for  stratified  bottoms, 
subroutines  which  are  being  developed  at  this  time. 

Several  subroutines  compute  coordinates,  travel  time  and  distance 
of  rays.  Calcomp  ray  plots  can  be  drawn.  The  intensity  is 
computed  such  that  each  ray  is  given  a  starting  intensity  according 
tc  the  directional  diagram  of  the  source.  This  intensity  is  being 
reduced  stepwise  by  medium  and  bottom  bounces.  A  number  of  rays, 
of  the  order  of  2000,  are  being  traced  in  that  manner  and  their 
intensity  parts  are  accumulated  at  given  distances  and  depth 
classes  incoherently. 

After  having  decomposed  the  input  set  P"  in  an  analogous  manner 
as  P1  the  internal  structure  of  the  first  function  h^  was 
demonstrated.  The  second  function  h'2  was  defined  after 
describing  the  sets  5’  and  A"  averaging  the  intensities  to 
sound  level  distributions  with  depth  and  total  propagation  losses. 

The  results  of  the  four  series  which  differed  most  were  demonstrated 
by  slides. 

The  influence  of  the  internal  waves  present  at  condition  4  brought 
in  the  experiment  a  double  value  of  total  propagation  loss  in  dB 
compared  with  the  loss  at  the  stationary  profile,  whereas  for  the 
computed  propagation  loss  a  rise  of  10$  occurred.  The  stationary 
profile  was  of  nearly  constant  gradient  with  a  sound  velocity  of 
1487  m/s  at  the  surface  and  of  1606  m/s  at  the  bottom.  Probably 
the  still  unknown  directional  diagram  of  the  source  and  reflection 
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behaviour  of  the  bottom  in  the  experiment  is  responsible  for 
this  disagreement. .  These  propertied  will  be  measured  soon 

and  put  into  the  program. 

With  this  simple  example  the  formal  concept  was  J«t 
But  a  concept  of  that  kind  will  really  be  helpful  x»  explain 
handling  and  controlling  the  problem,  when  large  numbers  of 
oceanographic  parameters  are  to  be  taken  into  account. 


APPLICATION  OF  RAY  TRACING  WITH  HORIZONTAL  GRADIENT 
TO  MONOSTATIC  BOUNDARY  REVERBERATION 


by 

L.B.  Palmer 

Naval  Research  Laboratory 
Washington,  D.C.,  U.S. 


ABSTRACT 

Presented  is  the  current  work  being  done  at  the  Naval  Research 
Laboratory  (NRL)  on  the  development  of  a  series  of  computer 
programs  to  predict  long  range,  low  frequency,  monostatic  boundary 
reverberation.  The  emphasis  is  on  the  ray  tracing  technique  and  its 
application  to  the  special  problems  of  estimating  the  transmission 
loss  of  rays  that  hit  the  boundaries  of  the  medium.  The  ray 
tracing  technique  is  to  increment  a  ray  from  point  to  point  along 
its  ray  path  by  evaluating  Taylor  series  expansions  in  arc  length 
of  various  ray  parameters  such  as  range,  depth,  travel  time, 
and  range  angle,  which  are  based  on  the  ray  equation.  Possible 
horizontal  variations  in  sound  speed  are  accounted  for  by  allowing 
multiple  input  profiles.  Also,  a  linearly  segmented  ocean  bottom 
and  a  flat  surface  are  assumed.  Monostatic  boundary  reverberation 
is  estimated  by  means  of  a  range  dependent  formulation  developed 
at  NRL.  An  underlying  assumption  of  this  formulation  is  that  when  a 
ray  encounters  a  boundary,  it  continues  to  propagate  in  the 
direction  of  specular  reflection,  while  a  small  amount  of  the 
incident  radiation  is  scattered  in  all  directions.  By  reciprocity, 
scattered  energy  will  return  to  the  sour ce-receiver  back  along 
those  ray  paths  emitting  from  the  source  and  passing  through  the 
boundary  point  at  which  the  hit  occurred. 
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INTRODUCTION 


The  purpose  of  this  talk  ^.s  to  present  the  current  work  being  done 
at  the  NRL  on  the  development  of  a  series  of  computer  programs  to 
predict  long  range,  low-frequency  monostatic  boundary  reverberation. 
The  emj  "\asi  b  will  be  on  the  ray  tracing  technique  and  its  application 
to  the  special  problems  of  estimating  the  transmission  loss  of 
rays  that  hit  the  boundaries  of  the  medium,  especially  the  bottom. 

The  tracing  of  rays  utilizes  an  iterative  technique  first  described 
in  Hudson  Laboratory  Report  150  by  W.A.  Hardy  et  al  in  1968. 

This  ray  tracing  model  accounts  for  the  possibility  of  a  horizontal 
sound  speed  gradient  by  allowing  for  multiple  input  profiles,  and 
assumes  a  linearly  segmented  ocean  bottom.  This  model  has  been 
programmed  in  FORTRAN  by  J.J.  Cornyn  of  NRL  as  part  of  a  series  of 
ray  tracing  and  transmission  lo&s  programs.  The  ray  tracing  model 
has  since  been  extracted  from  this  series  of  programs,  modified, 
and  then  adapted  to  the  purpose  of  calculating  monostatic  boundary 
reverberation.  Also,  computer  programs  have  been  written  which 
organize  the  ray  tracing  results  and  estimate  reverberation  by  means 
of  a  formulation  developed  by  J.T,  Warfield  of  NRL. 

I.  ORGANIZATION  OF  PROGRAMS 

The  operation  of  the  series  of  programs  proceeds  in  four  basic 
stages  [Fig.  1].  First,  since  the  modelling  of  the  velocity  field 
is  a  problem  unto  itself  because  of  the  possible  presence  of 
horizontal  variations  in  sound  speeds,  it  is  tackled  separately. 
Therefore,  first  a  magnetic  tape  is  generated  on  which  is  written 
the  multiple  profiles  located  at  discrete  ranges  along  the  track  of 
interest.  Each  profile  consists  of  the  sound  speed  and  its  first 
two  derivatives  with  respect  to  depth,  specified  at  discrete  depths 
and  the  ocean  bottom  as  discrete  depth  versus  range  points. 
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This  -tape  serves  as  input  to  the  ray  tracing  program.  This  program 
employs  the  iterative  ray  tracing  technique  to  determine  all 
boundary  hits  for  the  rays  of  interest.  These  points  as  well  as 
the  associated  travel  times,  ray  angles  at  the  boundary,  and 
transmission  losses  are  written  onto  an  output  tape.  An  initial 
ray  trace  may  not  furnish  enough  information  to  adequately  describe 
where  the  boundary  hits  occur.  If  this  is  the  case,  additional 
rays  will  be  traced  and  additional  ray  tracing  output  capes 
created. 

The  next  computer  program  uses  as  input  the  output  tapes  from  one 
or  more  ray  traces  and  reorganizes  this  information  onto  another 
magnetic  tape.  By  reorganizing  is  meant  the  creation  of  "order 
contours",  which  will  be  discussed  in  a  moment.  This  reorganizing 
procedure  is  done  separately,  due  to  computer  storage  requirements 
and  for  the  convenience  of  handling  multiple  ray  traces. 

Finally,  the  reorganized  tape  of  ray  tracing  information  is  the  input 
to  reverberation  calculating  program,  which  computes  the  boundary 
reverberation  function  at  discrete  times. 

II.  REVERBERATION  FORMULA 

Before  proceeding  to  a  description  >f  the  ray  tracing  model,  J.  would 
like  to  briefly  outline  the  model  used  for  predicting  monostatic 
boundary  reverberation. 

Equation  1  is  the  formula  used  for  determining  the  reverberation 
from  a  boundary  which  is  detected  at  the  source-receiver  at  time  t. 

R(t)=Ii  §  f  [  £  £  C(x)  /,  v  *  I  (p,  r,  x,  z)  ]  xdx  [Eq.  1] 

“0  pem(x)  rem(x)  pr''  * 

I)  is  the  source  level  relative  to  1  yd  and  the  remaining 
expression  represents  ai  tegration  over  the  three-dimensional 


boundary,  asstuning  separable  source-receiver  beam  patterns  and  an 
azimuthal  symmetry  of  the  boundary. 

An  underlying  assumption  in  the  development  of  the  reverberation 
formula,  is  that  when  a  ray  encounters  a  boundary,  it  continues  to 
propagate  in  the  direction  of  specular  reflection,  while  a  small 
amount  of  the  incident  radiation  is  scattered  in  all  directions. 

By  reciprocity,  scattered  energy  will  return  to  the  source-receiver 
back  along  those  ray  paths  emitting  from  fche  source  and  passing 
through  the  boundary  point  at  which  tl  .  hit  occurred. 

Consider  now,  the  expression  within  the  brackets,  which  is 
evaluated  at  each  range  x.  The  set  m(x)  is  a  set  indexing  all 
possible  ray  paths  from  the  source-receiver  which  pass  through  the 
boundary  at  range  x.  A  single  term  of  the  double  sum  represents 
the  contribution  to  reverberation  from  energy  travelling  from  the 
source  to  the  boundary  point  (x,  z)  along  path  p  and  returning 
along  path  r  [see  Figr  2], 

C  is  a  characteristic  function  which  is  equal  to  one  if  the 
reverberation  contribution  I,  attributed  to  the  ordered  ray  pair, 
(p,r),  will  be  detected  at  the  source-receiver  at  time  t,  and 
zero  otherwise.  The  times  of  detection  will  be  those  times  t 
3  vtisfying  the  inequality  of  Eq.  2. 

t-D  ^  Tp(x, z) +  Tr (x, z)  -  t  .  [Eq.  2] 

D  is  the  signal  duration  and  T  is  the  one-way  travel  time  for  a 
signal  travelling  tc  the  boundary  at  range  x,  along  the  path 
denoted  in  the  subscript. 

The  contribution,  if  detected,  is  given  by  Eq.  3. 

B(eo)  B(0_) 

I(P, -,*,*)  =  Lp(T’)Z)  ^  °C5p,\)  •  [Eq-  3] 
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B  and  B  are  the  transmitting  and  receiving  beam  patterns  respectively s 
and  the  L*s  are  the  transmission  losses  suffered  by  travelling  along 
the  respective  paths  and  o  is  the  boundary  scattering  strength. 

IIL  ORDER  CONTOURS 

In  order  to  get  a  visual  picture  as  to  how  the  possible  paths  between 
the  source-receiver  and  the  boundary  at  range  x  are  found,  an 
explanation  of  order  contours  is  necessary.  For  a  fixed  ray  R 
such  as  shown  in  Fig.  3,  each  boundary  hit  and  turning  point 
(both  of  which  are  called  occurrences)  is  assigned  a  positive 
integer  order.  Surface  hits  and  down  turning  points  (called  crests) 
are  defined  to  be  of  odd  order,  while  bottom  hits  and  up  turning 
points  (called  valleys)  are  of  even  order.  The  first  bottom  hit, 
or  valley,  for  the  fixed  ray  R  is  assigned  an  order  of  2. 

Subsequent  bottom  hits  and  valleys  are  assigned  4,  6,  8  and  so  on. 
Surface  hits  and  crests  are  assigned  orders  of  1,  3,  5,  etc.,  in 
such  a  manner  that  the  orders  of  the  occurrences  of  the  ray  R 
increase  with  range. 

Now  consider  the  contours  of  curves  determined  by  the  boundary 
hits  of  identical  order,  plotted -on  an  initial  source  angle  versus 
range  coordinate  system,  as  shown  in  Fig.  4*  Incidentally,  a 
given  contour  need  not  necessarily  be  a  continuous  curve,  but 
rather  several  continuous  and  disjoint  segments. 

For  the  reverberation  from  a  given  boundary,  only  those  contours 
corresponding  to  that  boundary  are  pertinent.  For  the  purpose  of 
illustration,  let  us  assume  that  bottom  reverberation  is  being  sought. 
To  determine  all  ray  paths  between  the  source-receiver  and  the 
bottom  at  range  x,  consider  a  horizontal  slice  at  range  x,  as 
shown  in  Fig.  5*  Denote  the  three  paths  having  initial  angles 
cf  0i,  83,  03  by  1,  2,  3  respectively.  Then  the  set  m(x)  consists 
of  the  integers  1,  2,  and  3,  and  hence  there  are  33,  or  9  possible 
routes  from  the  source-receiver  to  the  bottom  at  range  x  and 
back  again. 


In  reality,  the  computer  ray  tracing  program  will  produce  only  a 
finite  number  of  bottom  hits,  or,  in  other  words,  a  finite  number 
of  points  on  the  order  contours.  Figure  6  illustrates  two  order 
contours  of  the  type  produced  by  the  computer  program.  Here  the 
BT's  represent  the  bottom  hits  and  the  V’s  the  valleys  determined 
by  tracing  rays  with  initial  source  angles  of  8a,  0^,  6c,  etc. 

The  t;rue  order  contours  are  approximated  by  linearly  connecting 
the  known  bottom  hits  of  the  same  order.  Travel  time,  transmission 
loss,  and  ray  hit  angle  values  between  computed  hits  are  found  in 
interpolating  linearly  with  respect  to  range.  The  valley  type 
turning  points  are  displayed  to  indicate  the  end  points  of  the 
contours . 

IV.  RAY  TRACING  PROGRAM 

The  purpose  of  the  ray  tracing  then,  is  to  determine  the  boundary 
hits  and  the  transmission  loss,  travel  time,  and  ray  hit  angle  for 
e.-'.ch  hit.  The  rays  traced  should  be  those  rays  which  hit  the 
boundary  and  are  not  masked  out  by  the  source  beam  pattern.  Also, 
the  resulting  boundary  hits  should  approximate  the  true  order 
contours  as  accurately  as  possible.  In  the  search  fo^  a  ray 
tracing  model,  it  was  felt  that  such  a  model  should  allow  for 
multiple  nontrivial  sound  speed  profiles  and  ai.  irregular  bottom 
contour.  Such  a  model  was  developed  at  the  Hudson  Laboratories  at 
Columbia  University  of  New  York.  This  ray  tracing  technique  is  an 
iterative  one  in  which  a  ray  is  incremented  from  point  tc  point 
along  its  ray  path.  This  is  accomplished  by  means  of  evaluating 
Taylor  series  expansions  in  arc  length  of  various  ray  parameters 
such  as  range,  depth,  travel  time,  and  ray  angle.  These  expansions 
are  derived  from  the  ray  equation  shown  as  Eq.  4. 


A  £v(x,z7‘  “  grad  ' 


[Eq.  4] 


The  sound  speed,  v(x, z),  is  assumed  to  be  known  at  every  range 
and  depth  z,  throughout  the  two-dimensional  medium,  and  ds  is 
a  differential  increment  along  the  ray  path,  ?. 

The  present  version  of  this  ray  tracing  model  has  the  following 
features: 

1.  Accommodates  multiple  profiles  for  possible 
horizontal  variations  in  sound  speed. 

2.  Each  profile  is  defined  at  discrete  depths  with 
weighted  parabolic  interpolation  used  between 
specified  depths. 

3.  Assumes  a  linearly  segmented  ocean  bottom  and  flat 
surface,  with  specular  reflection  applied  to  both. 

4.  Uses  incremental  finite  Taylor  series  approximations 
to  ray  paths. 

5.  Transmission  loss,  ray  angle,  and  travel  time  computed 
at  boundary  hits  as  opposed  to  predetermined  ranges. 

6.  Each  ray  assigned  two  subliminal  companion  rays  for 
computing  the  amount  of  geometric  spreading. 

7.  Multiple  bottom  loss  tables  are  applied  to  various  range 
intervals . 

8  *  Individual  rays  are  terminated  upon  exceeding  input 

naximum  allowable  travel  time  or  bottom  loss  ceilings. 

9.  Up  to  500  rays  may  be  traced,  not  including  companion 
rays. 

Operationally,  rays  are  selected  whose  initial  source  angles  are 
from  that  part  of  the  source  beam  pattern  which  will  experience 
bottom  hits.  To  each  of  these  "primary"  rays,  is  assigned  two 
"companion"  rays  whose  initial  source  angles  tightly  bracket  that 
of  the  primary  ray.  The  rays  are  then  traced  through  the  medium 
which  is  described  by  the  multiple  profiles,  linearly  segmented 
bottom  and  flat  surface.  They  are  traced  on  an  individual  basis 


between  predetermined  "rectification"  ranges.  After  all  the  rays  are 
traced  to  one  of  these  ranges,  which  necessarily  include  the  ranges 
at  which  the  profiles  are  specified,  the  accumulated  boundary  hit 
and  turning  point  information  is  written  onto  a  magnetic  tape. 

A  computer  storage  area  is  then  reinitialized  and  the  ray  tracing 
continued  to  the  next  rectification  range  or  until  some  predetermined 
maximum  range  is  rea  :hed.  However,  individual  primary  rays  will 
be  terminated  prematurely  if  they  exceed  input  travel  time  or 
bottom  loss  ceilings.  For  each  boundary  hit  by  a  primary  ray,  six 
statistics  are  recorded. 


1.  Initial  angle  of  ray. 

2.  Order  of  hit. 

3.  Range  of  hit. 

4.  Travel  time. 

5.  Ray  angle  relative  to  the  boundary. 

6.  Transmission  loss. 


For  the  problem  to  which  the  present  model  is  applied.,  transmission 
loss  is  assumed  to  be  made  up  entirely  of  bottom  loss  and  geometric 
spreading  loss.  Bottom  loss  is  input  to  the  ray  tracing  program  as 
several  bottom  los j  versus  incident  angle  tables  which  are  applied 
to  different  range  intervals,  and  geometric  spreading  loss  is  assumed 
to  be  given  by  Eq,  5. 


1  1 

ascos  80  d80 

3 

a  d(sin  90)  ^ 

j  r  = 

x  cos  8x  dz 

x  cos  Qx  dz 

! 

j  where 

| 

1 

j  a  = 

unit  distance 

[Eq.  51 


range  of  boundary  hit 
depth  of  boundary  hit 
initial  source  angle  of  ray 
ray  angle  at  boundary  hit 
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The  companion  rays  [sec  Fig.  7]  are  used  to  calculate  the  factor  dz, 
by  determining  their  respective  depths  and  directions  at  the  range 
where  the  primary  ray  hits  the  boundary.  The  companion  ray  which 
has  not  yet  hit  the  boundary  in  this  region  and  is  still  being 
directed  towards  it,  is  used  bo  determine  dz.  After  a  primary 
ray  has  been  specularly  reflected  off  the  botbom,  the  amount  of 
boctom  loss  will  be  added  to  the  ray’s  accumulated  bottom  loss. 

V.  RAY  TRACING  TECHNIQUE 

Now  let  us  look  at  the  ray  tracing  technique  itself.  When  a  ray’s 
location,  as  well  as  its  ray  angle  and  travel  time  to  that  point, 
are  known,  the  sound  speed  and  its  spatial  derivatives  at  that 
position  are  used  to  determine  the  ray’s  location  and  associated 
parameters  at  a  point  further  along  the  ray  path.  A  single  profile 
could  be  employed  throughout  the  track,  but  often,  actually 
encountered  velocity  fields  require  the  specification  of 
different  profiles  along  the  track  of  interest.  Each  sound  speed 
profile  is  specified  at  discrete  depths  and  the  first  two  derivatives 
of  sound  speed  with  respect  to  depth  are  approximated  at  the  given 
depths  by  means  of  weighted  difference  equations.  Weighted 
parabolic  interpolation  is  applied  vertically  and  linear 
interpolation  horizontally  to  make  sound  speed  and  its  first  two 
depth  derivatives  functions  of  both  range  and  depth.  Therefore, 
the  first  derivative  of  sound  speed  with  respect  to  range  can  be 
approximated  at  any  range  and  depth  by  a  linear  first  difference  with 
respect  to  range.  With  the  use  of  this  horizontal  gradient  in 
the  Taylor  series  expansions,  it  is  hoped  to  reckon  with  the  effects 
on  ray  paths  due  to  horizontal  variations  in  sound  speed. 

Once  the  sound  speed,  its  first  two  derivatives  with  respect  to 
depth,  and  its  first  derivative  with  respect  to  range  are 
determined  for  a  known  location  on  the  ray  path,  various  convergence 
tests  are  performed  to  determine  how  far  to  increment  the  ray  path 
location,  and  associated  parameters,  by  evaluating  the  Taylor 
series  expansions.  Normally,  an  input  incremental  step  size,  A  , 


will  be  used  in  evaluating  the  Taylor  series  expansions.  However, 
if  certain  convergence  tests  fail,  the  step  size  will  be  f 

appropriately  reduced.  One  of  these  tests  require  that  the  ray  i 
path  does  not  experience  more  than  some  predetermined  maximum 
amount  of  bending  between  computed  ray  locations.  By  "bending" 
is  meant  the  change  in  the  trigonometric  sine  of  the  ray  angle. 
Another  test  requires  that  the  chosen  step  size  be  used  to 
predict  the  sound  speed  at  the  new  ray  location  within  a 
predetermined  accuracy  of  the  "true"  sound  speed  at  that  point. 

By  the  "true"  sound  speed  at  a  point  is  meant  that  value  found 
by  interpolating  between  input  profiles.  Regardless  of  these 
tests,  the  step  size  is  never  allowed  to  become  smaller  than  some 
predetermined  minimum  allovf.ble  step  size.  This  is  done  to  expedite 
the  incrementing  of  the  ray  path.  Therefore,  first  the  initial 
maximum  step  size,  A  ,  is  used  to  determine  a  trial  sine, 
given  in  Eq.  6. 


sin  0t  =  sin  80  +  a(Amax,  x0,  z Of  Z0,  D0,  G o  )  •  [Eq.  6] 


The  sound  speed  derivatives  Z,  D  and  G, 

Z  =  Bv/bz 
D  =  a3v/az2 
G  =  Bv/dx 


are  all  evaluated  at  the  known  ray  path  location,  (x0, zq).  A  new 
step  size,  A* ,  is  then  determined  from  Eq,  7. 


r 


A*  =  min  \ 


S  v. 


cos  0O  z0  6 


A 

max 


A 

max 


[Eq.  7] 
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where  S  is  an  input  maximum  allowable  change  in  sine,  say  0.02, 
Next,  using  the  new  step  A’,  a  new  location  (x^.,  z^.) ,  is  founu 
further  along  the  ray  path  by  evaluating  their  Taylor  series 
expansions. 

By  interpolation  between  two  known  profiles,  a  sound  speed 
is  determined  for  this  location.  Then  a  computed  sound  speed,  v  , 
is  found  for  this  location  using  the  Taylor  series  expansion  of 
Eq.  8. 

Dq  s 

vc(xt,zb)  =  v0(x0,z0)  +  Ze(zt  -  z0}  +-y  (zt-z0)  +G0(xt-x0)  [Eq.  8] 

If  the  computed  and  interpolation  sound  speeds  differ  by  more  than 
an  input  accuracy  criterion,  e,  then  a  new  step  size  is  found  from 
Eq;  9. 


A”  -mins  I^t-V°^ 


[Eq.  9] 


The  step  size  will  be  further  reduced  to  a  predetermined  minimum 

allowable  step  size  A  .  if  the  cosine  of  the  predicted  ray 

mxn 

angle  0^.  is  greater  than  one,  where  the  cosine  is  found  from 
Eq.  10. 

cos  0t  =  ctvt  [Eq.  10] 

and  where  b  is  found  by  evaluating  a  Taylor  series  in  arc 
length . 


The  final  step  size  is  never  allowed  to  be  less  than  A .  .  ,  except 

mm 

in  such  cases  as  approaching  a  boundary  or  a  known  profile 
location.  Once  the  final  step  size  is  determined,  a  new  ray  path 
location  and  corresponding  ray  parameters  are  found  by  evaluating 
the  different  series  for  this  step  size.  The  process  is  then 
repeated  with  the  new  location  being  taken  as  the  known  ray  path 
location. 
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Although  A  and  A  .  are  inputs  to  the  program,  they  may  be 
insx  mm 

redefined  as  a  ray  is  traced.  If  an  increase  results,  A  .  will 

mxn 

be  taken  to  be  the  vertical  distance  between  the  known  ray  depth 
and  the  nearest  depth,  in  the  direction  the  ray  is  headed,  at  which 
a  sound  speed  is  specified  in  the  input  profiles.  A  ^  will  be 
reduced  appropriately,  if  a  straight  line  projection  of  the  ray 
direction,  of  length  A  ,  passes  outside  the  medium.  Also,  when 
a  ray  location  is  determined  which  is  within  about  100  m  of  a 
boundary,  the  input  A  ^  will  be  used  as  the  iterative  step  size, 
without  the  various  tests  being  conducted. 

The  boundaries,  themselves,  are  assigned  tolerance  gates  of  some 
vertical  distance,  such  as  a  half  a  meter,  within  which  a  ray  is 
considered  as  having  struck  that  boundary.  In  some  cases, 
interpolation  will  be  employed  between  computed  ray  path  locations 
within  and  .^without  the  medium  in  order  to  determine  the  step 
size  necessary  to  increment  the  ray  to  the  boundary. 

The  described  ray  tracing  technique  represents  an  adaptation  of  the 
original  version  to  the  computation  of  transmission  loss  at 
boundary  hits  as  opposed  to  predetermined  receiver  ranges. 

CONCLUSIONS 

After  the  reverberation  and  ray  tracing  models  were  programmed, 
two  test  cases  were  run  to  check  the  validity  of  the  programs.  Both 
test  cases  involved  the  use  of  a  single  profile,  one  case  being  a 
constant  sound  speed  profile,  the  other  a  constant  gradient  profile. 
The  computer  results  for  these  two  test  cases  compared  quite  well 
with  the  corresponding  analytic  solutions.  The  ray  tracing  technique 
itself,  has  been  verified  for  a  multitude  of  cases. 

Next,  the  programs  were  run  on  a  case  for  which  measured  reverberatio 
data  existed.  This  case  entailed  a  track  of  fairly  long  range  over 
a  relatively  rough,  upward  sloping  bottom.  The  low-frequency  signal 
duration  was  about  30  s,  with  the  source  beam  pattern  dominated 
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by  the  main  lobe.  Input  boundary  scattering  strength  tables  were 
based  on  data  appearing  in  the  literature,  as  were  the  bottom  depth 
and  bottom  loss  tables. 

A  comparison  of  the  predicted  bottom  reverberation  versus  the  total 
measured  reverberation  is  shown  in  Fig.  8.  For  this  low  frequency 
study,  surface  reverberation  is  assumed  to  be  doppler  shifted. 

The  vertical,  or  reverberation,  axis  has  a  scaling  of  10  dB  between 
the  indicated  tick  marks.  The  relative  level  of  the  two  curves 
do  not  agree  as  well  as  had  been  hoped,  but  the  appearance  of  the' 
peaks  at  the  recorded  times  seem  to  correspond  fairly  well.  The 
peaks  of  the  predicted  curve  correspond  to  the  main  beam  of  the 
source  striking  the  bottom.  However,  the  relative  level  of  these 
peaks  are  largely  determined  by  a  small  angular  span  of  rays  whose 
initial  source  angles  usually  vary  by  less  than  1°.  These  so-called 
"crucial"  rays  are  those  rays  which  experience  both  valley  type 
turning  points  and  bottom  hits  with  very  low  grazing  angles  as  they 
propagate  down  range.  Some  of  these  rays  will  be  completely  RSR 
until  an  initial  bottom  hit  occurs  far  down  range,  at  which  point 
the  intensity  of  the  given  ray  will  dominate  those  of  other  rays 
hitting  the  bottom  nearby  which  will  have  accumulated  bottom  loss 
by  that  range.  These  "crucial"  rays  will  also  tend  to  encounter 
the  bottom  at  very  low  grazing  angles  where  the  bottom  loss  tables 
are  the  most  questionable.  Another  problem  is  the  occasional 
focusing  of  the  "crucial"  rays  near  the  bottom.  That  is,  they  are 
sometimes  in  near-caustic  regions  as  they  approach  the  ocean  bottom. 
Because  of  the  so-called  "crucial"  rays,  the  predicted  bottom 
reverberation  curve  is  sensitive  to  slight  changes  in  bottom  depth 
definition,  as  was  discovered  when  different  approximations  to 
the  ocean  bottom  were  used. 

The  linear  segmentation  of  the  bottom  is  also  a  problem  area.  The 
abrupt  changes  in  slope  between  linear  segments  tend  to  artificially 
break  up  the  wavefronts.  This  shows  up  as  radically  jagged  sections 
in  the  order  contours,  which  tends  to  weaken  one's  belief  in  the 
computer  approximation  of  the  contours  and  the  validity  of  the 


interpolation  performed  between  computed  bottom  hits.  Finally, 
averaging  over  reverberation  estimates  for  several  different 
headings  did  not  significantly  improve  the  results. 

Despite  these  apparent  problem  areas  the  main  difficulty  seems  to 
be  in  predicting  the  level  of  the  reverberation  curves.  Work  is 
continuing  on  this  problem  and  the  extension  to  the  bi static  case 
Investigation  is  underway  at  NRL  to  approximate  ray  intensities 
in  near  caustic  regions  and  to  more  efficiently  model  the  ocean 
bottom.  Also,  the  r~y  tracing  program  is  being  modified  to 
improve  the  iterative  procedure  and  the  spreading  loss  estimation 

DISCUSSION 

The  author  said  that  the  rough-surface  reflection  coefficient  was 
an  empirical  function  of  the  angles  of  incidence  and  reflection. 
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ABSTRACT 

The  purpose  of  this  paper  is  to  consider  the  efforts  performed  in 
the  development  of  ray  acoustic  propagation  analyses  from  the  view¬ 
point  of  advanced  sonar  system  design.  In  particular,  those 
parameters  required  of  such  analyses  for  the  enhancement  of  sonar 
system  design  and  performance,  are  identified  and  elaborated  upon. 

A  simple  ray  acoustic  shallow  water  propagation  model  is  used  as  a 
vehicle  for  identifying  and  relating  to  a  number  of  areas  in  which 
propagation  results  may  be  utilized  to  significantly  improve 
operational  sonar  system  effectiveness.  In  the  development  of  these 
examples,  the  spatial  distribution  of  the  acoustic  energy  together 
with  fluctuations  in  the  received  signal  energy  receive  prime 
consideration.  Expanding  upon  this  simple  basis,  a  sonar  design 
concept  is  developed  by  analytically  incorporating  propagation 
features  that  can  be  provided  by  suitable  models. 


The,  orientation  of  my  discussion  today  is  somewhat  different  from 
those  that  have  been  presented  so  far  in  this  conference.  I  intend 
to  talk  not  as  an  underwater  acoustician  actively  involved  in 
propagation  analysis  and  the  development  of  ray  trace  programs. 
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Rather,  from  my  view  point  as  an  engineer  faced  with  the  problem  of 
developing  sonar  systems  to  operate  with  maximum  efficiency  in  the 
ocean  medium,  I  wish  to  present  some  of  my  requirements  for 
information  on  propagation  phenomena. 

My  reasons  for  attempting  such  a  presentation  are  to  outline  to 
representatives  of  the  propagation  community  those  aspects  of 
propagation  which  bear  directly  on  design  decisions,  and  to 
indicate  why  simple  measures  of  the  energy  lost  by  a  signal 
propagating  through  the  sea  are  not  in  themselves  sufficient 
inputs  to  the  sonar  design  process.  In  doing  this,  I  hope  to 
indicate  to  you  some  specific  design-related  parameters  which 
characterize  propagation  phenomena  and  about  which  relatively 
little  information  exists.  Specifically,  it  is  my  objective  to 
illustrate  why  the  spatial  and  temporal  behaviour  imposed  on  signals 
by  the  medium  present  inputs  to  the  sonar  design  process  equal  in 
importance  to  a  knowledge  of  the  propagation  loss  between  two  points. 

I  hope  that  this  overview  will  encourage  studies  in  which  attention 
is  focused  on  developing  more  information  about  the  mechanisms 
governing  these  spatial  and  temporal  characteristics  which- currently 
are  often  "washed  out"  or  ignored  in  analyses  directed  at  developing 
numbers  to  characterize  average  propagation  loss. 

The  general  points  I  wish  to  consider  are  outlined  in  Fig.  1.  The 
usual  goals  of  a  programme  to  develop  propagation  models,  using  either 
ray  tracing  techniques  or  normal  mode  analysis,  are  twofold, 
namely: 

1.  To  delineate,  on  the  basis  of  velocity  profiles  and 
other  environmental  data,  the  dominant  paths  via  which 
acoustic  signals  propagate  from  source  to  receiver. 

2.  To  develop  a  means  for  estimating  the  average  acoustic 
power  lost  by  a  signal  as  it  propagates  between  two  points. 
The  loss  mechanisms  involve  accounting  for  geometrical 
spreading  of  the  wavefront,  ref lection/scattering,  and 
absorption  effects  within  the  medium  or  at  its  boundaries. 


Once  these  goals  are  accomplished,  the  sonar  designer  can  readily 
estimate ; 

1.  The  average  value  of  the  signal  present  at  his  array;  and 

2.  The  depths  at  which  he  can  expect  high  and  low  levels 
of  signal  power. 

In  general,  however,  a  comprehensive  sonar  design  process  requires 
more  detailed  knowledge  concerning  the  influences  of  the  medium 
upon  the  signal.  Specifically,  the  designer  would  like  to  have 
at  his  finger  tips: 

1.  The  spatial  distribution  of  the  signal  arrivals  at 
the  array . 

2.  The  statistics  of  both  spatial  and  temporal  fluctuations 
imposed  on  the  signal  by  the  medium;  and 

3.  In  the  case  of  short  term  (transient  or  pulse)  signals, 
the  "impulse  response"  of  the  medium;  i.e.,  a  measure  of 
the  time  and  frequency  dispersion  imposed  on  the 
original  signal . 

For  the  remainder  of  this  paper,  I  intend  to  focus  on  the  three 
areas  listed  above,  and  illustrate,  by  some  simple  examples,  how 
knowledge  of  these  phenomena  may  be  incorporated  into  the  sonar 
design  process.  In  particular,  I  wish  to  concentrate  on  the  active 
sonar  problem  since  this  area  is  most  urgently  in  need  of  more 
sophisticated  techniques. 

As  a  starting  point,  I  would  like  bo  consider  a  simple  example  to 
illustrate  how  a  knowledge  of  the  spatial  distribution  of  signal 
arrivals  can  impact  on  the  design  of  a  sonar  array.  The  specific 
case  of  interest  here  is  one  where  the  impact  is  greatest,  namely, 
a  vortical  array  for  use  in  a  shallow  water  environment.  As  used 
here,  the  term  shallow  water  implies  those  instances  where  propagation 
takes  place  over  ranges  very  large  relative  to  water  depth,  and 
where  the  water  depth  itself  is  large  compared  with  the  wavelength 
of  the  signal  being  propagated. 


In  such  instances  (particularly  when  ranges  are  such  that  little 
energy  arrives  via  a  direct  path),  the  signal  arrives  via  multiple 
surface  and  bottom  reflections.  These  multiple  arrivals  imply 
that  the  dominant  signal  energy  is  distributed  over  a  range  of 
arrival  angles  rather  than  being  concentrated  at  a  single  angle. 

If  this  range  of  angles  is  known,  then  the  array  can  be  designed 
to  spatially  "match"  this  angular  window,  and  thus  insure  that 
maximum  signal  power  is  captured  by  the  array. 

This  can  easily  be  illustrated  by  use  of  a  very  simple  ray  trace 
model,  which  describes  the  gross  features  of  propagation  in  a 
shallow  water  environment.  This  model  is  similar  to  the  one 
devised  by  McPherson  and  Daintit-h  [Ref.  1],  and  elaborated  on 
by  Smith  [Ref.  2].  The  salient  features  of  this  model  (which  is 
based  on  averaging  over  ray  cycles)  are  shown  in  Fig.  2, 

The  details  of  the  model  on  which  the  results  to  be  presented  are 
based  can  be  found  in  the  references.  Basically,  it  treats  the 
large  number  of  rays  involved  in  transporting  energy  from  source 
to  receiver  as  a  statistical  ensemble.  The  major  features  of  the 
propagation  mechanisms  are  preserved,  but  by  use  of  spatial  and 
temporal  averaging  over  the  ranges  and  depths  involved,  fine 
details  (i.e,,  Lloyd’s  mirror  interference  phenomena)  are  smoothed 
out.  However,  the  fundamental  postulates  and  assumptions  which 
guide  the  analysis  are  indicated  in  the  figure. 

The  postulates  are  self  explanatory;  the  assumptions  are  found  to 
be  reasonable  withir.  the  objectives  of  the  analysis,  namely,  a  gross 
prediction  of  the  propagation  loss  which  preserves  the  basic 
sensitivity  of  the  result  to  the  dominant  environmental  parameters. 
These  are; 

1.  Water  depth- 

2V  Sound  speed  profile. 

3.  Surface  loss. 

4.  Bottom  loss. 
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Although  sound  velocity  profiles  are  seldom  observed  to  have  equal 
velocity  at  all  depths,  the  isovelocity  model  often  predicts 
results  close  to  those  observed  experimentally  in  shallow  water. 
The  most  significant  departure  from  isovelocity  behaviour  occurs 
for  strong  surface  ducts  and  severe  upward  or  downward  refraction. 
The  model  can,  however,  be  readily  extended  to  include  these  cases 

The  geometry  for  a  simple  isovelocity  situation  is  shown  in 
Fig.  3,  along  with  the  results  of  the  analysis. 

In  the  isovelocity  case  £  ,  the  loss  per  cycle  from  boundary 
reflection  for  a  ray  with  initial  angle  0  ,  is  found  to  be 

£  =  bg  u  +  bb  0  =  b(e)  (sin0  =  B,  cos  0=1) 


where 


b  = 


The  results  of  Smith’s  analysis  [Ref.  2]  show  that  the  ray  cycle 
method  gives  the  following  value  of  transmission  ratio  r(r) 
under  the  isovelocity  assumption.  (The  transmission  ratio  is 
the  ratio  of  the  intensity  at  range  r  to  that  at  1  yd)  and  is 
given  by: 


T(r) 
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The  expression  in  the  figure  contains  twice  the  integral  from  0 
(horizontal)  to  ..  limiting  angle  9f,  For  the  symmetrical 
situation,  thi  is  equivalent  to  including  all  rays  from  **0£  to 
+0£»  Note  that  this  integral  is  of  the  form  of  the  area  under  a 
Gaussian  (normal)  curve..  In  their  paper,  McPherson  and  Daintith 
[Ref,.  1]  obtained  a  similar  j  sovelo-'  -ty  result  by  analysing  the 
number  of  br-unces,  rather  than  working  with  ray  cycles.  They 
showed  that  the  result  is  the  same  for  the  sloping  bottom  case  if 
average  depth  is  used  for  D,  In  addition,  they  derived  equations 
that  are  quite  similar  to  those  shown  here  for  propagation  loss 
under  negative  gradient  conditions. 


The  distribution  of  signal  arrival  angles  as  determined  from  the 
equation  for  the  transmission  ratio  per  radian  is  shown  in  Fig.  4« 

Note  that  in  isovelocity  water  the  received  energy  is  normally 
distributed  in  arrival  angle  about  a  mean  angle  9  =  0  with  a 
standard  deviation  c  that  increases  with  increasing  depth  and 
decreases  with  increased  range  and  bottom  absorptivity.  The  upper 
limit  of  the  integral  for  an  omnidirectional  receiver  is  tt/2 
corresponding  to  the  two  cut-off  points  shown  in  the  figure. 

The  implications  on  array  design  are  immediately  evident.  First, 
if  one  designs  a  vertical  array  to  have  a  beamwidth  of  6g  between 
say,  the  6  dB  down  points,  whan  steered  to  the  maximum  response 
axis,  the  choice  of  9g  should  be  dictated  by  the  distribution  of 
signal  arrivals  as  shown  in  the  figure.  These,  in  turn,  depend 
upon  the  water  depth,  the  range  between  source  and  receiver,  and 
the  combined  bottom  and  surface  losses.  While  the  first  two 
parameters  are  usually  definable  for  the  conditions  over  which  an 
array  must  operate,  the  determination  of  the  surface  and  bottom 
losses  is  not  so  easily  made.  However,  it  is  interesting  bo  note 
that  sufficient  information  for  the  designer  is  obtained  by 
developing  estimates  of  average  bottom  and  surface  losses  over  the 
ocean  area  of  interest.  Thus,  a  useful  effort  would  be  to  continue 
work  aiajed  at  developing  estimates  of  average  bottom  and  surface 
losses,  parameterized  in  terms  of  controlling  physical  mechanisms 
(i.e.,  windspeed,  etc.)  in  areas  of  interest  to  the  sonar  community. 


The  preceding  example  results  in  conclusions  which  may  seem  somewhat 
obvious  to  people  versed  in  propagation  analysis.  However,  this 
exmple  example  xv as  intended  mainly  to  illustrate  the  coupling  which 
exists  between  the  design  process  and  a  knowledge  of  medium  effects} 
the  two  cannot  be  separated,  nor  can  simple  estimates  of  energy  loss 
suffice.  In  this  case,  the  spatial  distribution  of  the  energy 
arrivals  impact  on  design  choices}  in  turn,  the  design  impacts  on 
guiding  the  propagation  analysis  by  delineating  the  environmental 
parameters  of  importance. 
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This  analysis  also  points  out  the  types  of  shallow  water  propagation 
models  of  use  to  the  sonar  designer.  Essentially,  he  looks  for  a 
model  which  provides  maximum  visibility  of  the  physical  parameters 
of  the  ocean  environment  which  influence  sonar  performance. 

In  this  context,  a  simple  approximate  model,  which  can  be  rapidly 
executed  may  be  favoured  over  an  elaborate  ray  trace  model  which 
requires  large  computer  facilities.  At  the  same  time,  the  designer 
is  aware  of  the  limitations  of  such  simple  approximate  models. 

In  the  one  employed  here,  for  example,  the  use  of  incoherent 
energy  addition  of  ray  paths  precludes  accounting  for  spatial 
phase  cancellation  effects  among  rays  (e.g.,  the  Lloyd's  mirror 
effect).  However,  as  long  as  the  designer  is  aware  of  the  model 
limitations,  he  can  use  various  sub-models  to  investigate  "fine 
structure"  effects  as  required.  Basically,  this  amounts  to  an 
iterative  application  of  propagation  models;  i.e.,  begin  with  the 
simplest  physical  model  available,  and  then  expand  the  analysis 
through  the  use  of  sub-models  as  the  design  development  dictates. 

Of  even  more  relevance  to  the  sonar  designer  are  data  concerning 
the  characteristics  of  fluctuations  imposed  on  the  signal  by  the 
medium.  Sound  fluctuations  in  the  ocean  are  observed  for  virtually 
every  type  of  propagation;  rface  ducts,  deep  sound  channels,  sea 
surface  or  bottom  reflected  paths  —  to  name  a  few.  All  of  these 
propagation  paths  can  be  characterized  by  an  observed  mean  value  of 
propagation  loss  around  which  instantaneous  values  of  transmitted 
energy  are  distribute 

Too  often,  modelling  efforts  are  focused  on  predicting  the  mean 
values,  with  little  or  no  attention  given  to  estimation  of  the 
fluctuating  components .  Yet,  often  the  fluctuations  are  of  primary 
interest  to  the  sonar  designer,  in  that  they  impact  on  his 
selection  of  processing  times  and  his  estimations  of  the  spatial 
stability  of  the  signal  over  tho  array  aperture.  Recent  work,  such 
as  that  of  Nichols  and  Young  [Ref.  3],  and  Dyer  [Ref.  4]  have 
provided  considerable  information  on  the  character  and  statistics  of 
fluctuations;  however,  more  extensive  studies  are  needed  if  an 
adequate  data  base  is  to  be  established. 
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As  an  example  of  how  such  information  may  be  incorporated  into  the 
sonar  design  process,  consider  the  problem  of  estimating  the 
spatial  phase  stability  of  a  signal  over  the  aperture  of  a  sonar 
system.  Here,  it  is  important  to  understand  that  the  degree  to 
which  a  sonar  array  can  improve  signal-to-noise  gain  and  provide 
directional  information  is  directly  related  to  the  degree  of 
stationarity  of  the  phase  of  the  signal  wavefront  over  the  array 
elements.  When  the  array  is  "steered"  to  the  direction  from  which 
the  signal  arrives,  any  "jitter"  or  fluctuations  in  phase  among 
elements  of  the  array  can  cause  a  degradation  in  performance,  in 
the  sense  that  the  signal  is  not  perfectly  "in  phase"  at  every 
element  of  the  array. 

This  can  be  illustrated  by  considering  the  response  of  a  simple 
two-element  array  to  a  signal  arriving  from  an  angle,  8,  as  shown 
in  Fig.  5.  The  two  identical  receivers  each  generate  an  output 
voltage  e0  in  response  to  the  incoming  wave.  However,  the 
difference  in  signal  path  length  between  the  two  receivers  results 
in  a  spatial  phase  difference  between  their  outputs.  Referenced  to 
the  geometric  center  of  the  two-element  array,  the  sum  of  the 
voltages  from  the  two  hydrophones  is  given  as 


eT  =  2e0cos  sin  8)  =  2e0  cos  (u)T) 


where  k  =  ~  , 

A 


element  spacing. 


Since  k  =  u)/c,  where  uj  is  the  angular  frequency  and  c  the  velocity 
of  sound  in  the  medium,  the  phase  term  can  be  written  as 

,  aid  .  .  v 

cos(~  sm  8)  =  cos  U)T 

d  sin  8 

where  T  =  -  represents  the  time  delay  between  constant  phase 

c 

(spatial)  arrivals  at  each  hydrophone.  Note  that  when  the  signal 
arrives  at  an  angle  8  =  0  (i.e.  broadside),  there  is  zero  time  delay, 

and  the  output  voltage  of  the  array  is  a  maximum;  i.e.  2eo .  Suppose, 
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however,  that  fluctuations  in  the  medium  produce  a  fluctuation  in 
the  delay  time,  such  that  T  =  T0+  T(t)j  i.e,  the  delay  is 
characterized  by  a  mean  value  and  a  fluctuation  about  that  mean. 

In  this  case,  for  a  signal  incident  along  maximum  response 
axis  (8  =  0),  T0  =  0,  but  t(t)  may  be  finite.  The  output  of 
the  array  in  this  case  is  given  as 

eT  =  2e0  cos[(JJT(t)] 

where  T(t)  may  be  treated  as  a  random  variable. 

To  develop  an  average  response  we  are  interested  in  the  expected 
value  of  cos[o)r(t)].  This  value  is  given  as 

(cos[u>T(t)]  >  =  J  p(t)  COSUlTdT 

wher-e  p(t)  is  the  probability  distribution  of  the  time  delay 
fluctuation.  If  we  assume  that  the  values  of  T(t)  are  Gaussianly 
distributed,  and  have  a  zero  mean,  then 

1  -T2/ 202 

(cos[u)T(t)])  = -  cos  iut  e  dr 

V  2ttot  ^ 

The  result  of  the  integration  yielding 

-m2a2/2 

(cos[«JT(t)])  =  e 

Since  we  assumed  the  amplitudes  of  the  arrivals  to  be  constant, 
the  expected  value  of  the  array  output  voltage  is 

-tt)So2/ 2 

(© j )  2e0  e  . 

Note  that  as  the  array  response  is  maximum.  However,  as 

the  magnitude  of  the  fluctuations  increase,  as  reflected  in 
increasing  values  of  c®  of  the  value  of  (©^,)  decreases  from 
2e0  for  "on-beam"  arrivals. 


Intuitively,  one  might  expect  the  magnitude  of  tv  phase  fluctuations 

to  increase  in  proportion  to  the  separation  between  elements.  Thus, 

one.  might  assume  a  ad,  i.e.  a  =  yd,  where  v  is  the  constant 

T  T 

of  proportionality.  In  this  case,  the  expected  output  voltage 
of  the  ideal  two-element  array  is  given  as 

u)2d 

-y  —j— 

<eT>  =  2ee  e 

Note  that  the  degradation  in  response  increase  as  the  frequency, 
as  well  as  with  increased  element  spacing. 

The  result  derived  here,  although  for  a  simple,  if  not  trivial, 
case  is  illustrative  of  how  such  fluctuation  information  can  be  used 
by  the  sonar  designer.  The  analysis  techniques  can  be  readily 
extended  both  to  include  multi-element  large  arrays,  and  also  to 
allow  for  fluctuations  in  the  amplitude  of  signal  arrivals.  The 
main  point  here,  however,  is  the  indication  of  a  need  for  better 
measurements  of  characteristic  delay  time  fluctuations  induced  by 
the  medium,  and  most  importantly  the  development  of  propagation 
models  that  incorporate  and  predict  these  parameters. 

Up  to  this  point,  1  have  been  stressing  the  parameters  required  of 
propagation  models,  in  addition  to  signal  attenuation,  to  produce 
a  driving  influence  on  the  sonar  design  and  advanced  development 
process.  All  too  often  it  seems  that  the  signal  processing  people 
proceed  independently  of  the  environmental  people  with  the  result 
that  the  idealizations  that  appear  to  perform  so  well  in  radar 
fail  to  reach  expectations  in  the  underwater  medium.  My  message 
then  is  that  although  the  question  of  the  best  mathematical  fit 
to  a  sound  velocity  profile  is  important,  the  character  of  the 
outputs  in  both  space  and  time  that  are  obtained  from  a  propagation 
model  are  of  perhaps  greater  importance. 
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As  a  final  example  to  illustrate  this  point,  I  wish  to  consider  a 
problem  in  the  active  sonar  area.  In  general,  as  the  threat  becomes 
more  quiet,  active  sonar  must  greatly  increase  its  capability  to 
enable  more  likely  detections  to  occur  in  a  shorter  time  frame j 
thereby  reducing  the  threat  of  counter  detection  and  localization. 

To  achieve  this  increase  in  capability  in  the  face  of  reverberation 
it  appears  that  more  than  just  an  increase  in  source  level  will  be 
required.  In  particular,  it  seems  reasonable  that  increased  gains 
in  signal  processing  must  be  obtained.  These  gains  in  turn  can 
only  result  from  increased  knowledge  of  the  environmental  conditions 
and  their  effect  on  the  transmitted  signal. 

Proceeding  from  this  point,  let  us  conceptualize  a  hypothetical 
system  that  might  be  realizable  in  the  not  too  distant  future. 

Such  a  system  is  shown  in  block  diagram  form  in  Fig.  6.  We 
hypothesize  an  advanced  active  sonar  sub-system  capable  of 
transmitting  and  receiving  fanar  signals  and  displaying  the  processed 
results.  In  addition,  we  include  an  environmental  measurement 
sub-system  capable  of  providing  real  time  environmental  parameters 
including  measurements  of  the  sound  velocity  profile.  Interfacing 
these  two  sub-systems,  we  envision  a  ray  trace  propagation  sub-system 
providing  on-site  estimates  of  the  signal  propagation  situation  based 
on  the  environmental  (and  historical)  inputs.  The  outputs  of  the 
propagation  sub-system  perhaps  with  an  operator  interface  dictate 
the  type  of  sonar  transmission  mode  to  be  employed,  e.g.  bottom  bounce 
direct  path,  etc.,  for  the  mission  to  be  performed  and  also  the  type 
of  signal  to  be  transmitted.  In  addition,  the  propagation  sub-system 
also  determines  the  type  of  receiver  processing  both  spatial  and 
temporal  to  be  employed. 

To  illustrate  how  these  sub-systems  might  interact,  let  me  consider 
another  simple  example,  based  on  some  of  the  propagation  modelling 
that  was  alluded  to  earlier  in  the  talk.  We  assume  an  environment 
as  depicted  in  Fig.  7  which  might  be  representative  of  some 
generalized  area  during  the  winter /early  spring  months  where  upward 
refraction  is  obtained  just  below  the  surface  because  of  the  lack 


of  surface  heating.  Under  the  assumption  of  isogradient  conditions, 
one  characteristic  of  the  propagation  behaviour  is  that  of  a  deep 
surface  duct  as  shown  in  the  figure.  If  a  SVP  gradient  of  0.025/s 
is  assumed  then  at  short  ranges  the  limiting  ray  of  the  duct  would 
reach  a  vertex  of  about  200  ft.  The  rays  "trapped"  within  this 
limiting  ray  would  be  surface-reflected  and  their  energy 
contribution  would  arrive  after  that  of  the  limiting  ray.  The 
transmission  angle  of  the  limiting  ray  would  be  much  le.  s  than  6°, 
implying  that  the  surface  reflected  rays  are  all  reflected  at 
grazing  angles  less  than  this  value.  Since  surface  reflected 
paths  will  undergo  very  little  scattering,  at  reasonable  sea 
states  the  sea  surface  reflection  may  be  considered  essentially 
specular.  Under  these  conditions,  the  surf ace- reflected  rays  will 
all  arrive  at  essentially  the  same  intensity  as  that  of  the 
non-ref lected  limiting  ray. 

The  travel  time  for  the  limiting  ray  is  the  smallest  among  the 
duct  arrivals  and  consequently  the  energy  travelling  this  path 
is  the  first  to  arrive  from  the  duct.  The  energy  travelling  e 
remaining  ducted  paths  will  arrive  in  a  non-uniform  manner  following 
that  of  the  energy  travelling  the  limiting  or  refracted-only  ray. 

The  order  of  arrival  follows  the  number  of  surface  reflections. 

In  other  words,  the  energy  travelling  the  one  surface  reflection  path 
will  be  following  by  that  travelling  the  two  surface  reflection  paths 
and  so  on.  In  terms  of  our  system  design,  it  is  envisioned  that 
this  type  of  information  would  be  automatically  displayed  by  our 
hypothesized  ray  path  computer. 

As  an  approximation  for  illustrating  the  method  by  which  the  receiver 
processing  might  be  determined  by  the  ray  path  computer,  let  us 
assume  that  the  ducted  multipaths  arrive  uniformly  in  time  with  an 
equal  delay  of  &  seconds  between  each  path.  If  the  arrivals  are 
assumed  of  equal  intensity  the  normalized  receded  signal  can  be 
expressed  as  shown  by  the  equation: 


s(t)  =  rr  £  cos(wb  -  &i)  . 
N  i=l 


In  this  equation  the  effects  of  doppler  scattering  causing  shifts 
in  the  received  frequency  are  neglected.  The  energy  envelope  of 
the  received  signal  can  easily  be  expressed  by  noting  that  it  is 
of  the  form  of  the  off-axis  response  of  an  equally  spaced  coherent 
array  [Ref.  5] •  In  particular,  the  envelope  of  s(t)  is  given  by 

S3  =  rsin^ui4^2)  3 
e  lN  sin  (  U)A/  2fJ 

In  general  there  will  be  a  large  number  of  arrivals  from  the  duct* 
and  the  delay  between  arrivals  can  be  considered  to  approach  zero. 
This  fact  may  be  incorporated  in  the  expression  by  letting  N*»«* 
while  hr* 0  with  NA^T,  where  t  is  the  total  time  delay  of  the 
duct.  Under  these  conditions 


«c  - 


N-»«,  A**0,  NA"*T  . 


This  expression  is  termed  the  coherency  factor  and  represents  the 
normalized  energy  arriving  via  the  duct  as  a  function  of  the  transmitted 
frequency  and  the  duct  tine  delay. 


At  this  point,  we  require  some  knowledge  of  the  environment 
parameters.  From  the  ray  trace  sub-system,  we  are  able  to  determine 
a  value  for  the  average  duct  time  delay.  If  we  use  the  shallow 
water  propagation  model  referred  to  earlier  in  the  talk,  this  is 
indeed  possible.  In  fact,  Smith  has  shewn  that  for  a  shallow  water 
isogradient  duct,  the  average  duct  time  delay  can  be  expressed 
approximately  by 

T  -  i  (|)  (-if-5)  . 

In  this  equation  R  is  the  range  of  interest,  cg  and  c^  are 
the  sound  velocity  values  at  the  source  and  the  vertex  of  the 
limiting  ray  respectively,  and  c  is  the  nominal  velocity  of  sound 
in  sea  water  (5000  ft/s).  For  a  SVP  gradient  of  0.025/s,  a  depth 


difference  of  150  ft  between  the  source  and  the  vertex  of  the 
limiting  ray  and  for  the  purpose  of  illustration  a  nominal  range 
of  6.5  kyd,  then  a  time  delay  of  1  ms  results. 

Given  this  computation  by  the  propagation  sub-system,  the  sonar 
operator  may  determine  that  the  energy  transmitted  in  the  duct  at 
least  at  short  ranges  is  quite  concentrated  in  time.  Thus  if  his 
sonar  can  transmit  pulse  burst  transmissions  which  would  generally 
have  pulses  long  compared  to  1  ms  he  need  not  concern  himself  with 
time  dispersion  caused  by  the  surface  duct  transmission.  In  other 
words,  in  this  situation  he  concludes  that  he  need  not  worry  about 
the  time  dispersion. 

In  addition,  if  all  of  the  reverberation  energy  associated  with  the 
termination  of  a  given  pulse  decays  before  the  arrival  of  the  next 
transmitted  pulse,  then  the  pulse  period  within  a  burst,  i.e.,  the 
pulse-on  to  pulse-off  duty  cycle,  should  be  reasonably  stable.  All 
of  these  results  arise  from  the  fact  that  for  a  1  ms  duct  delay 
as  much  as  a  100$  perturbation  caused  by  medium  variations  amounts 
to  only  a  maximum  2  ms  d«J  y  in  the  signal  time. 

At  this  point,  the  operator  has  used  the  propagation  results  to 
indicate  the  type  of  propagation  to  be  encountered  and  the  extent 
of  the  time  delay  and  in  turn  has  used  this  information  to  determine 
the  type  of  signals  he  should  transmit.  Because  of  the  short  time 
delay,  he  has  concluded  that  gated  pulse  signals  will  not  be 
dispersed  in  time  enough  to  cause  them  to  be  distorted  greatly. 
Frequency  distortion  caused  by  the  multipath  situation  is  another 
matter  of  interest  to  the  operator.  If  he  has  a  range  of  transmission 
frequencies  available  to  him  he  would  like  to  choose  a  value  to 
minimize  the  fading  characteristics  of  the  signal. 

Since  the  coherency  factor  has  previously  been  determined  this 
function  with  the  appropriate  time  delay  obtained  from  the  propagation 
results  may  be  used  by  the  operator  to  choose  from  among  the 
frequencies  available  to  him. 
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Thus,  based  on  calculations  provided  by  the  propagation  computations 
the  operator  may  select  within  the  available  range,  the  transmission 
frequency  for  which  the  signal  fading  is  minimized. 

The  total  propagation  situation  in  our  example  consists  of  a 
bottom  bounce  mode  in  addition  to  the  surface  duct  transmission. 

[See  Fig.  7].  In  general,  the  energy  arriving  via  these  two 
paths  will  interfere  and  may  cause  degradation  in  the  total 
received  energy.  However,  with  an  adequate  ray  trace  program  in 
the  system,  it  is  conceivable  that  an  on-the-spot  determination  of 
the  arrival  time  difference  between  the  two  paths  may  be  made  and 
that  this  information  could  be  used  to  either  select  an  appropriate 
processing  mode  to  eliminate  the  effects  or  by  properly  assessing 
the  characteristics  to  use  them  to  advantage.  In  the  former 
approach,  one  might  select  a  processor  similar  to  the  Reke  System 
in  radar  which  attempts  to  separate  the  various  energy  arrivals  and 
then  to  coherently  recombine  them  to  enhance  the  total  received 
energy.  In  the  latter  technique,  by  properly  assessing  the 
propagation  situation,  the  operator  might  be  able  to  range  on  the 
basis  of  received  signal  strength. 

Let  me  consider  this  latter  technique  in  more  detail  since  it 
provides  a  good  illustration  of  incorporating  propagation  information 
into  the  operational  problem.  Without  exact  knowledge  of  bottom 
topography,  it  is  always  difficult  to  compute  the  bottom  bounce  paths; 
however,  a  reasonable  approach  is  to  assume  an  essentially  flat 
bottom  over  the  bounce  region.  If  this  assumption  is  made,  then  for 
the  isogradient  case  representative  calculations  indicate  that  the 
time  delay  of  the  single  bottom  bounce  path  relative  to  the 
limiting  ray  (for  a  gradient  of  0.025/s)  may  be  computed  to  be  in  the 
range  of  10  ms.  This  result,  of  course,  is  extremely  sensitive  to 
parameter  changes.  However,  a  good  lower  bound  may  bi  8  ms  with 
an  upper  bound  as  high  as  50  ms. 

Incorporating  this  result  with  the  surface  channel  delay  of  1  ms 
indicates  that  the  total  impulse  response  of  the  measurement 


channel  in  the  short  range  case  we  are  considering  consists  of  the 
narrow  surface  duct  response  plus  a  bottom  arrival  impulse  (probably 
a  series  of  bottom  arrivals)  delayed  8  ms  to  SO  ms  with  respect  to 
the  surface  arrivals.  To  keep  this  example  simple  let  us  neglect 
the  multiple  bottom  bounce  arrivals  on  the  basis  of  energy 
contributions. 

In  the  duct,  the  incoherent  energy  is  given  by  the  incoherent 

addition  of  the  energy  of  each  multipath  arrival?  E^.  This  implies 
that  is  given  by  the  number  of  significant  arrival  paths  times 

the  energy  per  path. 


EIC  =  NEp  * 

The  coherent  energy,  on  the  other  hand,  is  given  by 

E_  =  N2  E  a  . 

L  pc 

In  this  equation  a  is  a  measure  of  the  coherency.  If  a  is  1, 

c  c 

the  energy  arrivals  are  coherent  and  maximum  energy  is  obtained. 
Eliminating  N  from  the  two  equations  results  in  the  equation  shown 
for  the  coherent  energy. 

E 

Ci  r?  '  • 

C  E 

P 

If  one  bottom  arrival  is  considered,  then  it  contributes  the  energy 
per  path  E  modified  by  the  bottom  loss  factor  p.  Thus,  the 

X** 

ratio  of  the  bottom  to  duct  energy  is  given  by 


RBD  = 


e /e 


IC 


aJEr 


Since  fche  source  is  common,  the  RBD  ratio  can  be  expressed  in 
decibels  in  terms  of  the  duct  and  bottom  path  transmission  losses. 
This  expression  is  given  below.  It  expresses  the  ratio  of  the  bottom 
to  duct  energy  in  terms  of  the  duct  and  bottom  transmission  losses, 
the  duct  coherency  factor  and  the  bottom  loss  parameter 


RBD  =  2(11^  -  TLd)  -  10  lg  dn  -  10  lg  p 
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To  utilize  this  expression,  we  determine  duct  and  bottom  transmission 
losses  from  the  shallow  water  (range  much  greater  than  depth) 
propagation  loss  model  of  Smith  [Ref.  2].  An  approximate  expression 
for  the  transmission  loss  in  an  isogradient  duct  is 


e~aR<1\ 

TLd“101  ' 


In  this  equation  Rd  is  the  horizontal  range  between  the  source 
and  the  receiver,  ft  is  the  volumetric  attenuation  coefficient  as  a 
function  of  frequency,  L  is  the  depth  of  the  duct  measured  from 
the  surface  to  the  limiting  ray,  and  8^  is  the  angle  of  the 
limiting  ray  at  the  source  which  is  determined  by  the  equation 

2(cb-c0)  l/a 

0*  -  i-r; — a  • 

Using  a  sound  speed  gradient  of  0.025/s  in  these  equations  and  an 
assumed  duct  depth  of  200  ft,  we  determine  that 


e, 
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In  addition,  for  the  200  ft  duct  we  determine  the  duct  transmission 
loss  as 

TL^  =  -28  -  10  lg  Rd-  4.34ftRd  • 

For  the  bottom  transmission  loss  assume  spherical  spreading  with  a 
range  of  essentially  the  same  as  that  of  the  duct  and  a  volumetric 
attenuation  q  .  Thus,  the  simple  spherical  spreading  loss  is 


TL 


b  =  -20  lgRd-  4.34ftRd 


From  our  previous  calculations  a  representative  value  of  the 
coherency  factor  may  be  determined  to  be  -28  d3.  Thus,  we  may 
combine  our  results  to  obtain  the  expression  for  the  ratio  of  the 
bottom  and  duct  energy. 


RBD 


2(28-  10  lg  R ,)  +  28  -  10  lg  p  . 
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Choosing  a  short  horizontal  range  of  6 . 5  kyd, 

RBD  =  6  -  10  lg  0  . 

A  reasonable  bottom  loss  value  in  the  range  4  dB  to  9  dB  may  be 
chosen  fr;m  Urick  [Ref.  6]  for  a  sand  type  composition.  Using  this 
range  of  values  yields  the  range  of  RBD  values  of  from  0  db  to  -3  dB. 
Thus,  on  the  basis  of  these  calculations  it  may  be  expected  that 
the  duct  and  bottom  energy  arrivals  are  of  the  same  order  of 
intensity.  This  is  a  significant  result  for  the  sonar  operator  in 
that  it  implies  that  pulse  spreading  caused  by  the  bottom  arrival  is 
of  the  same  order  of  intensity  as  that  of  the  duct  arrivals. 

The  combination  of  the  duct  and  bottom  arrivals  at  a  receiver  implies 
a  time  spreading  of  an  individual  pulse  and  the  generation  of  an 
interference  pattern  in  the  region  of  time  overlap.  The  magnitude 
of  this  interference  is  most  simplyg/investigated  by  the  simple 
Lloyd's  Mirror  type  of  pattern  generated  by  two  arrivals  of  the  same 
frequency.  A  plot  of  the  expected  envelope  in  the  region  of  overlap 
is  shown  in  Fig.  8  as  a  function  of  operating  frequency  and  the 
difference  between  the  oath  arrival  times.  The  curve  is  also 
parameterized  by  the  relative  intensities  of  the  two  arrivals. 

An  interpretation  of  this  diagram  is  as  follows.  Assuming  that  the 
energy  envelope  of  the  initial  pulse  arrival  is  at  zero  dB  (for 
convenience)  when  th^  energy  arrives  on  the  second  path,  an 
interference  pattern  is  established  during  the  duration  of  overlap  of 
the  two  pulses.  The  magnitude  of  the  envelope  during  this  overlap 
period  is  a  function  of  the  degree  to  which  the  arrivals  coherently 
add  (measured  by  fi-)  and  their  relative  intensities  (measured  in  dB) 
and  may  be  found  from  the  diagram.  When  time  progresses  beyond  the 
overlap  period,  the  envelope  of  the  response  reduces  to  that  of  the 
second  arrival.  In  terms  of  the  results  we  have  found  for  the 
surface  duct  and  bottom  path,  since  the  two  arrivals  are  at 
essentially  the  same  energy,  one  would  expect  a  relatively  continuous 
pulse  with  a  discontinuity  (either  an  increase  or  a  decrease 
depending  on  the  f t  product)  during  the  period  of  overlap.  From 


previous  calculations  the  difference  between  the  duct  and  bottom 
arrivals  can  be  expected  to  be  on  the  order  of  8  ms  to  50  ms. 

In  these  ranges  one  would  expect  both  increases  and  decreases 
during  the  overlap  of  the  duct  and  bottom  arrivals.  In  addition, 
of  course,  the  pulse  will  be  spread  in  time  by  the  later  bottom 
arrival . 

Although  the  computations  we  have  presented  may  seem  somewhat 
complicated,  in  reality,  of  course,  they  could  easily  be  programmed 
on  a  mini -computer.  In  this  manner,  the  operator  could,  at  least 
ideally,  be  presented  with  an  indication  of  the  type  of  pulse 
distortion  that  he  might  expect  and  could  either  alter  his 
frequency  or  signal  shape  accordingly. 

In  summary,  I  have  tried  to  indicate  the  importance  of  propagation 
models  from  the  view-point  of  sonar  design  and  future  sonar 
systems.  In  doing  so,  I  have  tried  to  emphasize  that  simple 
predictions  of  transmission  loss  are  not  sufficient  for  these 
purposes.  Rather,  more  emphasis  must  be  placed  upon  the  temporal 
and  spatial  influences  cf  propagation,  the  fluctuations  encountered, 
and  the  interference  effects  that  may  occur  in  transmitted  signals. 

I  hope  that  these  remarks  may  serve  to  stimulate  the  direction  of 
future  work  and  operation. 
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DISCUSSION 


■When  asked  whether  sonars  should  not  be  designed  fcr  more  general  use 
than  just  shallow-water  operation,  the  author  said  that  the  trend  was  toward 
more  specialised  sonars,  otherwise  the  necessary  compromises  become 
too  difficult. 


PROPAGATION  ANALYSIS:  USUAL 

OBJECTIVES 

X.  DELINEATE  DOMINANT  PATHS  OF  ENERGY  TRANSMISSION 

2.  ESTIMATE  LOSS  IN  ENERGY  AS  SIGNAL  PROPAGATION  FROM 
SOURCE  TO  RECEIVER 


PROPAGATION  ANALYSIS :  ADDITIONAL 

OBJECTIVES 

1.  SPATIAL  DISTRIBUTION  OF  SIGNAL  ARRIVALS  (MULTIPATH  EFFECTS) 

2.  STATISTICS  OF  SPATIAL  AND  TEMPORAL  FLUCTUATIONS  IMPOSED 
ON  SIGNAL  BY  MEDIUM 

3.  IMPULSE  RESPONSE  OF  MEDIUM  AS  IT  INFLUENCES  SHORT  TERM 
(PULSE  OR  TRANSIENT  SIGNALS) 


FIG.  1 


SIMPLE  SHALLOW  WATER 
PROPAGATION  MODEL 

REFERENCES:  SMITH.  P.W.  Jr.."SOUND  TRANSMISSION  IN  SHALLOW  WATER: 

PART  I:  ANALYSIS. "  BB  t  N  REPORT  NO.  1563.  OCT. 24. 1967 
McPherson.  j.d..  and  daintith.  mj..  "practical  model  of 

SHALLOW  WATER  ACOUSTIC  PROPAGATION."  J.ACOUST.  SOC.  AM.  VOL.  41 


FIG.  2 


RAY  CYCLE  MODEL  (AFTER  SMITH) 


POSTULATES  1.  ENERGY  TRAVELS  ALCNG  RAY  PATHS 

2.  ENERGY  IS  BOTH  CONTINUOUSLY  ATTENUATED  (VOLUME)  AND 
INTERMITTENTLY  ATTENUATED  (SURFACE  +  BOTTOM  REFLECTION/ 


SCATTERING) 

3.  ENERGY  TRAVELS  AT  LOCAL  SOUND  SPEED  ALONG  RAY 


ASSUMPTIONS  1.  ONLY  RAYS  IN  WATER  CONTRIBUTE  (NO  BOTTOM  REFRACTED  PATH) 

2.  ENERGY  SCATTERED  TOO  FAR  OUT  OF  BEAM  ESSENTIALLY  LOST 

3.  BULK  OF  ENERGY  RECEIVED  CLUSTERED  ABOUT  IDEAL  RAY  PATHS 

4.  SOURCE SEQULN  .  OF  IMPULSES.  RECEIVED  SIGNAL  p2  (tk  l.t. 
SUPERPOSITION  OF  ENEROY  ARRIVALS  (NO  INTERFERENCE  EFFECTS) 


LOSSES 

P  =  bs9  *  bb0 

b  =  bg+b(3 


ISOVELOCITY  CASE 

TRANSMISSION  LOSS 

b0  Nw  -  10  log  T 

where 
r 


bs «  surface  loss  per  ray  cycle 
per  radian 

bb=  bottom  loss  per  ray  cycle 
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SIGNAL  ARRIVAL  ANGLES 
FOR  ISOVELOCITY  CASE 


FIG.  4 
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WTERFERENCE  INTENSITY  RELATIVE  TO  WITIAL  ARRf/AL 
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FIG.  5 


However,  when  r  *  r0+  r(t) 
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°r*  ■  Variance  of  distribution  of  r(t) 


FIG.  6 


POSITIVE  ISOGRADIENT 
SHORT  RANGE  PROPAGATION 
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SURFACE  CHANNEL 


FIG.  7 


TWO  RATH  LLOYD  MIRROR  INTERFERENCE  PATTERN 


FIG.  8 
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DETERMINATION  OF  THE  INTENSITY  OF  SOUND  AT  ARBITRARY  POINTS 
IN  THE  SOUND-FIELD  OF  A  SOURCE  IN  A  HORIZONTALLY  LAYERED  MEDIUM 


M.J.  van  der  Scheur 
Physics  Laboratory,  T.N.O., 

The  Hague,  Netherlands 

The  following  method  is  based  on  finding  all  the  rays  through  a  point 
P  and  to  add  the  corresponding  intensities. 

Consider  a  layered  medium  with  a  linear  depth-velocity  profile.  ‘The 
ray-path  and  the  corresponding  intensity  losses  are  calculated  according 
to  a  number  of  well-known  formulae  (Fig.  1). 

To  determine  the  intensity  at  P,  we  add  an  extra  layer  through  P  to 
the  velocity/ depth  profile.  When  two  rays,  leaving  the  source  with 
starting-angles  close  to  each  other,  intersect  the  level  through  P  on 
both  sides  of  P,  there  will  be  generally  at  least  1  ray  with  a  starting 
angle  between  the  two  mentioned  ones,  which  reach  the  level  at  P, 

An  iterative  process  will  give  us  the  value  of  80 • 

Of  course  we  need  a  set  of  good  starting  values  for  this  iterative 

process.  Therefore,  we  define  the  characteristic  velocities  of  the 

depth- velocity  profile.  These  are  the  greatest  velocities  of  each 

layer,  provided  that  this  value  is  greater  than  the  velocity  at  the 

source-depth,  Cg,  and  greater  than  all  values  occurring  between  the 

source  and  the  layer  considered.  The  corresponding  characteristic 

Co 

values  of  K(80)  ~  •  gives  the  start  angles  of  the  rays,  which 

COS  Oq 

will  turn  back  just  at  the  limit  of  a  layer. 

The  horizontal  distance  between  the  source  and  a  point  on  the  ray,  S, 
is  a  function  of  K(80)  :  S  =  S(K(80)).  We  consider  one  specified 
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depth  level,  S  consists  of  a  number  of  pieces  of  the  type  A Sx,  ASa , 
AS3  (see  Fig.  1).  We  can  write  in  general 


Ax  Ag  A3 

s  =  A  <AS‘h  +A  usa).  +1s1  as3>. 

JQ 

Consider  the  derivatives  of  S  [Fig.  2],  We  see  that  can  be 

written  as  the  sum  of  two  monotonic  functions*  one  increasing  and  the 
other  decreasing.  From  this  we  can  conclude  that  in  the  interval 

jq 

between  two  characteristic  values  of  K(8),  ^  has  two  zero  points:  K(6*). 

Now  we  add  all  the  values  K(0*)  for  which  |j|?  =  0  to  the  array  of 
characteristic  values  K(0),  This  means  thau  1a  the  interval  between 
two  successive  values  of  array  K,  the  function  of  the  horizontal  distance 
S  at  a  certain  level  is  monotonic.  When  two  rays  with  successive 
values  K[i 3  andK[i+l]  intersect  the  considered  level  on  both  sides  of  P, 
there  will  be  exactly  1  ray  with  a  starting  angle  0  between  K[i]  and 
K[i+1],  intersecting  the  level  in  P. 

In  our  computer  model  the  ray  path  will  be  symmetric,  fixed  by  three 
values  a,  b,  c  [see  Fig.  3], 

As  you  can  see,  we  make  a  difference  between  the  direct  and  indirect 
”fcll 

rays.  The  n  intersection  with  the  level  through  P  of  the  ray  with 
value  K[i]  lies  at 


DISTANCE  (K[i],n)  »  a  +  [§]  *  b  +  [^]  *  c 


for  the  indirect  rays: 


INDIST  (K[i],n)  =  d  +  [|]  *  b  +  [2^]  *  c 

K[i]:  the  array  of  characteristic 
starting  values. 

For  every  value  of  n  we  decide  if  array  K[i]  must  be  completed  with 

JO 

values  for  which  ■35?  =  0. 

OK 


Define:  MAX(n)=  m*X  DISTANCE (K(i),n)  MIN(n)=  m^n  DlSTANCE(K(i),r?) 

When  for  the  horizontal  distance  between  the  source  and  the  point  P, 
range  P,  the  following  relation  holds: 

MIN(n)  <  range  P  <  MAX(n), 

then  there  will  be  rays  intersecting  the  level  considered  for  the  k 
time  in  P. 

When  we  repeat  this  process  for  n=i,  N  both  for  direct  as  well  as 

indirect  rays,  then  we  will  find  M  rays  going  through  P.  For  each 
of  these  rays  we  determine  the  intensity  and  finally  we  find  for  the 
total  transmission  loss  at  P 

Nsp,  -  1°  ^ 

Of  course  there  are  some  restrictions  in  the  present  computer  model. 
Except  for  the  restriction  of  a  linear  depth-velocity  profile  the  most 
important  assumption  is  that  the  sea  surface  is  a  flat  plate  in  order 
to  obtain  a  symmetric  ray  path.  However,  there  is  the  possibility 
of  giving  an  attenuation  factor  for  each  surface  reflection.  For  the 
bottom  similar  assumptions  are  made. 

When  we  introduce  a  waving  surface,  the  simplicity  of  the  computation 
disappears  since  the  function  S  becomes  much  more  complicated. 


DISCUSSION 

The  author  said  that  no  comparisons  had  yet  been  made  with 
measurements. 
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POSITION  AND  SHAPE  OF  THIS  SURFACE  SHADOW  ZONE 


B.  de  Raigniac 
SACLANT  ASW  Research  Centre 
la  Spezia,  Italy 


From  surface  reverberation  ^uaies  conducted  by  SACLANTCEN*s 
Target  Classification  Group*,  it  has  been  found  necessary  to  know 
before  any  sea  experiment  the  precise  position  of  the  first 
surface  shadow  zone.  Since  normal  ray  tracing  programs  usually 
do  not  answer  such  a  question,  a  simple  computer  program  has  been 
implemented  which,  on  a  small  shipborne  computer,  calculates 
within  seconds:  distance,  extent  and  maximum  thickness  of  the 
shadow  zone  as  a  function  of  the  source  depth. 

From  submarine  detection  trials  subsequently  performed,  it  was 
seen  that  a  more  accurate  shape  of  this  shadow  zone  was  desirable 
and  an  extension  of  the  previous  program  is  underway. 

DISTANCE.  EXTENT  AND  THICKNESS  OF  THE  SHADOW  ZONE 

Figure  1  shows  the  rays  which  limit  the  shadow  zone  both  in  range 
and  depth  when  the  source  is  between  the  minimum  velocity  depth  and 
the  critical  depth.  At  the  surface,  the  shadow  zone  is  bounded  by 
the  rays  which  have  zero  grazing  angle.  The  maximum  dupth  of  the 
shadow  zone  is  the  depth  at  which  a  ray  horizontal  at  source  becomes 
horizontal  again.  A  similar  figure  would  be  obtained  with  a  source 


W,  Bachmann  and  B.  dB  Raigniac,  "The  Calculation  of  the  Surface  Backscattering 
Coefficient  of  Underwater  Sound  from  Measured  Data",  SACLANTCEN  Technical  Memorandum 
No.  174,  November  1971,  NA7Q  UNCLASSIFIFD 
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between  surface  and  minimum  velocity  depth;  the  shadow  zone  is 
then  of  much  larger  dimensions.  Once  the  Snell's  constant  k  of 
the  ray  is  known,  the  distance  to  the  shadow  zone  is 


D  =  £  R±  (sinai_1  -  sir.  ,  [E<1*  *] 

where  R^  is  the  radius  of  ray  curvature  in  each  layer  i, 
ai~l  an£*  ai  are  Sraz^nS  angles  of  the  ray  at  the 
boundaries  of  the  ith  layer,  and  the  summation  extends  over 
s.11  layers  from  the  source  to  the  surface. 

The  extent  of  the  shadow  zone  is  given  by 

S  =  2  £  R.  (sin  a.  ,  -  sin  a.)  ,  [Eq,  2] 

X  1  l-A  1 

where  the  summation  is  now  taken  over  the  layers  between  the 
source  and  the  critical  depth  2^.  . 

The  maximum  thickness  of  the  shadow  zone  is 


=  Z 


j  V  “  V 


[Eq.  3] 


where  V  is  the  sound  speed  at  the  source  depth, 
s 


In  Fig,  2  these  three  quantities  are  plotted  as  a  function  of 
source  depth  fcr  a  typical  summer  Mediterranean  sound  velocity 
profile. 


SHAP F  OF  THE  SHADOW  ZONE 

Figure  3  gives  an  idea  of  the  shadow  zone  as  obtained  by  a 
conventional  ray  tx'acing  program.  A  precise  determination  of  the 
shape  requires  a  large  number  of  rays  and  therefore  a  large  amount 
of  computer  time  is  involved,  particularly  on  a  small  shipborne 
machine. 

An  alternative  method  is  to  look  directly  for  the  rays  (in  caustics 
or  limiting  rays)  which  delimit  the  shadow  zone.  The  exact  solution 
is  mathematically  difficult.  The  shape  is  then  approximated  as  the 
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locus  of  the  points  where  at  closely  spaced  depths  a  ray  vertexes. 

As  limiting  rays  may  occur  when  strong  negative  gradients  are 
predominant,  a  comparison  is  performed  at  each  depth  between  ranges 
of  the  vertex  points  and  of  intersection  of  rays  vertexing  at 
shallower  depth.  The  points  retained  are  those  which  indicate  a 
smaller  extent  of  the  shadow  zone,  as  indicated  in  Fig.  4. 

Figure  5  shows  hew  the  shape  of  the  shadow  zone  varies  as  a  function 
of  the  source  depth  for  the  same  summer  velocity  profile.  It  should 
be  noted  how  the  shape  can  be  complicated  and  how  the  single 
knowledge  of  distance,  extent  and  maximum  thickness  may  not  be 
sufficient . 

In  Fig.  6  different  shapes  are  obtained  from  a  typical  winter 
sound  velocity  profile.  It  should  be  noted  that  the  maximum 
shadow  zone  depth  given  by  Eq.  3  is  not  always  reached  because  of 
limiting  rays. 


DISCUSSION 


In  reply  to  a  question,  the  author  reiterated  the  point  that  a  complete 
ray  tracing  was  not  required  to  obtain  the  shape  of  the  shadow  zone} 
the  method  described  being  very  much  simpler  than  that. 


FIG.  1 
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INFORMAL  PRESENTATIONS 


1 •  OCEANIC-ACOUSTIC  EXPERIMENTS  AT  SACLANTCEN 
T.D.  Allan  (SACLANTCEN,  La  Spezia,  Italy) 

hast  year,  at  the  Director's  suggestion,  a  project  in  oceanic- 
acoustics  was  set  up  at  the  Centre  with  scientists  from  at  least 
3  groups  playing  some  part  in  the  planning  and  execution  of 
experiments.  The  project  represents  a  serious  effort  to  break 
down  the  older  concepts  of  ASW  oceanographers  and  acousticians 
following  broadly  the  same  research  lines  but  working  independently 
of  one  another  and  planning  separate  cruises.  With  oceanographers 
and  acousticians  working  together'  at  sea  the  acoustic  trials 
benefit  from  having  much  more  detailed  environmental  data  than 
usually  available  in  trials  while  any  prediction  of  sonar  conditions 
that  the  oceanographer  is  likely  to  make  from  the  environmental 
data  can  be  checked  by  the  acoustician. 

The  relevance  of  ray-tracing  to  this  type  of  research  is  obvious. 

It  provides  the  bridge  between  the  observed  sound  velocity  or 
temperature  structure  and  the  predicted  sonar  conditions  with  which 
to  compare  the  observed  sonar  conditions. 

The  first  oceanic-acoustic  cruise  was  carried  out  in  February  1970 
in  the  N.W.  Mediterranean  [Fig.  11. 

The  area  chosen  for  the  experiment  had  been  surveyed  the  previous 
February  during  a  multi-ship  expedition  sponsored  by  the  NATO  Sub¬ 
committee  on  Oceanographic  Research  (MEDOC  69),  and  it  appeared 
that  during  winter  the  Ligurian  Sea  could  be  divided  into  3 
different  zones  according  to  the  temperature  structures  shown  in 
Figs.  2-4.  Acoustic  runs  wore  therefore  planned  as  shown  by 
the  full  lines  in  Fig.  1.  Although  hampered  by  appalling  weather 
the  experiment  had  some  success.  It  certainly  showed  that 
acousticians  and  oceanographers  could  integrate  their  programmes 
with  little  or  no  mutual  interference. 


We  were  fortunate  to  be  offered  the  services  of  one  French  and 
two  Italian  minesweepers  besides  our  own  research  vessel  MARIA 
PAOLINA.  The  procedure  was  planned  as  follows: 

(a)  Over  the  three  shorten  profiles  (10  miles  each)  the 
receiving  ship,  MARIA  PAOLINA,  was  to  remain  stationary  at  the 
end  of  the  profile  with  hydrophones  placed  at  25m,  50m,  100 m, 

200m  and  500m  depth.  XBT 1  s  to  be  taken  every  10  minutes. 

(b)  The  French  minesweeper  to  act  as  firing  ship  opening 
and  closing  range  on  MPG,  firing  deep  and  shallow  (25m  and  200m) 

^  lb  explosive  charges  every  5  minutes. 

(c)  The  two  Italian  minesweepers  to  follow  behind  the 
firing  ship  at  distances  of  and  5  miles  taking  XBT ' s  every 
10  minutes  [see  Fig.  51. 

(d)  Each  run  to  be  repeated  3  times  in  succession. 

(e)  At  the  end  of  the  third  run  the  French  minesweeper 
to  transmit  CW  pulses  (50  ms  pulses  at  20  seconds  interval, 

3*5  kHz)  for  a  period  of  1  hour  at  a  fixed  range  of  10  miles. 

All  fou"  ships  to  continue  to  tr-ke  XBT's  every  10  minutes  during 
this  time-series  experiment. 

(f)  CW  transmissions  to  be  repeated  at  a  range  of  5  milts 
for  1  hour,  with  the  XBT  ships  again  aligned  between  source  and 
receiver. 

(g)  For  the  long  profile  (20  miles)  across  the  stable  and 
unstable  zones  the  same  procedure  to  be  used  but  with  XBT's  every 
5  minutes.  (It  might  be  noted  here  that  the  major  factor  in 
deciding  the  time  interval  between  XBT's  was  simply  the  total 
number  available  to  us  for  the  four  ships.) 

(h)  Two  self-recording  oceanographic  buoys  with  suspended 
arrays  of  thermistors  (20  from  the  surface  to  300  m)  to  be  situated 
at  positions  Bx  and  Bs  so  that  r-wo  further  temperature  inputs 
would  be  available  along  the  long  profile  2  and  one  each  along 
profiles  3a  and  3b. 
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Ae  stated,  the  weather  was  bad  f.n d  for  more  than  B0%  of  the  time 
the  minesweepers  had  to  lie  in  harbour.  Not  only  was  the  programme 
severely  curtailed  but  the  unusual  frequency  of  high  winds 
contributed  to  a  radically  different  oceanographic  environment 
in  197D  from  that  found  in  1969. 

Of  the  four  acoustic  profiles  planned,  only  one  complete  profile 
and  1/3  of  another  (the  latter  made  with  two  support  ships  instead 
of  three)  could  be  made.  This  was  a  great  pity  because  the 
procedure  went  very  smoothly  at  the  first  attempt.  The  locations 
of  the  XBT  measurements  are  shown  in  Figs.  6  and  7 • 

Data  Processing  and  Analysis 

All  XBT  traces  were  digitised  and  hence  interpolated  to  provide  n 
temperature  reading  every  2 m  down  to  300  m  and  every  20  m  from 
300  m  to  500  m,  Each  trace  was  preceded  by  an  identifying  ship 
number,  and  distance  along  the  profile  from  the  receiving  ship, 
MARIA  PA0LINA. 

Through  the  good  offices  of  Dr  J.B.  Hersey  and  Capt  P.  Wolff,  the 
Fleet  Numerical  Weather  Central,  Monterey  agreed  to  carry  out  the 
i ange-dependent  ray-tracing  required  to  compare  the  computed  and 
observed  rays  from  the  two  sources  to  the  five  hydrophones. 

The  following  information  was  requested,  firstly,  for  a  ray-trace 
made  with  only  the  BT  taken  from  the  firing  ship  and  subsequently 
for  a  ray-trace  made  with  4  simultaneous  BT 1 s  taken  along  the 
profile  at  the  moment  of  fire: 

For  each  ray: 

(i)  Path  length  in  metres, 

(ii)  Travel,  time  in  milliseconds  to  the  third  decimal  place, 
if  possible. 

(iii)  Propagation  loss  anomaly,  in  decibels  defined  as 
A  =  10  log10 

/V 

where  I  is  acoustic  intensity  at  receiver 


I0  is  acoustic  intensity  at  lm  from  the  source 
is  the  straight-1 Lne  distance  from  source  t  > 
receiver  in  metres 
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(iv)  Angle  at  the  source,  in  degrees  measured  to  the 
horizontal . 

(v)  Grazing  angle,  in  degrees,  of  the  surface-reflected 
rays  for  each  reflection. 

This  involved  a  considerable  number  of  ray  tracings  —  over  a 
hundred  for  the  one  profile  successfully  completed  as  planned. 

Although  at  the  tj.me  of  writing  the  data  analysis  is  well  underway 
it  is,  perhaps,  too  early  to  try  to  draw  conclusions.  There  have 
been  some  delays  caused  by  working  with  two  computers  of  different 
generations  but  for  most  of  the  events  we  now  have  plots,  for 
each  of  ten  hydrophone- source  combinations,  of  the  various  ray 
paths  predicted  by  four  (or  more,  by  moving  backwards  or  forward 
in  time)  temperature  traces  and  by  the  single  temperature  trace 
at  the  firing  ship.  Similarly  all  of  the  observed  acoustic  data 
are  analysed  and  plotted  for  each  event.  Within  a  few  months  we 
should  be  able  to  say  how  well  they  match. 

A  few  words  on  another  oceanic-acoustic  project  carried  out  again 
with  the  help  of  Monterey.  This  is  an  investigation  of  the  acoustic 
effect  of  an  oceanic  front. 

Such  a  front  was  surveyed  in  December  1970  by  0.  Jcha-.nessen  using 
the  thermistor  chain  on  the  LEE  (US) .  Temperature  traces  measured 
every  6  minutes  along  a  45-mile  track  thro  Ii  the  rront  were 
digitised  and  taken  to  Monterey  together  with  one  detailed 
section  of  8  miles  in  the  centre  of  the  frontal  region  where  the 
sampling  interval  was  90  seconds. 

From  these  traces,  propagation  loss  has  been  calculated  in  both 
directions  through  the  front  (east-west  and  west-east)  for  a  variety 
of  source  and  receiver  depths  and  for  three  frequencies.  The 
computer  calculation  already  show  that  there  are  significant 
differences  in  propagation  loss  in  the  two  directions. 

In  December  1971  it  is  hoped  to  measure  the  effect  directly  in 
another  joint  oceanic-acoustic  cruise. 


FIG.  1  LOCATION  OF  ACOUSTIC  RUNS  PLANNED  FOR  THE  FEBRUARY  CRUISE 
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FIG.  2  TEMPERATURE  PROFILES  TAKEN  APPROXIMATELY  EVERY  HOUR  AT  ONE  STATION  IN  THU 
STABILISED  COASTAL  ZONE  (McD0C  69) 
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FIG.  5  PLAN  FOR  THE  OCEANOGRAPHIC/ACOUSTIC  MEASUREMENTS  ALONG  THE  SHORT  RUNS. 

Each  ship  to  take  an  X8T  every  10  minutes  with  the  firing  ship  exploding  deep  ond  shollow  chcrges 
every  5  minutes.  MARIA  PAOLINA  stationary  at  the  end  of  the  profile  with  hydrophones  at  5  depths 
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ACOUSTIC  PROPAGATION  THROUGH  OCEANIC  FRONTS 
M.J.  Daintith  (A,U,W.  E.,  Portland,  Dorset,  U.K.) 


A  rdsumd  was  given  of  recent  experiments  in  the  Mediterranean, 
East  of  Malta,  designed  as  a  preliminary  test  of  the  effect  of 
oceanic  fronts  on  acoustic  propagation.  It  was  apparent  that 
in  certain  circumstances  the  propagation  loss  was  considerably 
increased  by  the  presence  of  a  front j  and  it  seems  likely  that 
the  front's  presence  could  lead  to  significant  bearing  errors. 


SUMMARY  AND  RECOMMENDATIONS  -  A  PERSONAL  VIEW 


by 

Richard  H.  Clarke,  Conference  Secretary 


An  unfortunate  first  impression  might  be  that  ray  tracing  is 
dull,  simply  because  the  idea  is  a  very  old  one.  In  fact,  ray 
tracing  is  the  single  most  powerful  method  available  to  us  for 
understanding  the  complexities  of  sound  propagation  in  the 
ocean.  And  gaining  this  understanding  is  quite  literally  all 
our  business,  I  will  take  the  main  topics  of  the  conference  in 
turn,  attempt  to  give  a  personal  summary  of  what  has  been 
presented,  and  suggest  where  future  research  should  be  directed. 

The  Ocean 

We  learned,  or  were  forcefully  reminded* that  fc r  from  being 
uniform  in  the  horizontal  direction  there  are  drastic  changes 
that  occur  across  fronts  (analogous  to  atmospheric  fronts) , 
and  that  these  fronts  can  be  found  in  all  oceans  and  seas. 
Further  than  that,  the  temperature  and  salinity  (and  hence 
the  sound  velocity)  are  not  continuous  in  depth  but  often 
proceed  in  steps  of  isothermal  and  isohaline  water,  separated 
by  thin  sheets  across  which  the  temperature  and  salinity 
change  abruptly. 

As  acousticians  we  must  thank  the  oceanographers  for  their 
awkward  findings,  but  press  them  to  give  us  more  detail;  both 
about  the  fine  structure  of  temperature  and  salintiy  and  its 
global  occurrence , 
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Sound  Velocity 


That  vital  link  which  joins  the  oceanographer* s  description  of 
the  ocean  to  the  acoustician's  has  made  great  strides  in  the 
last  thirty  years  or  so.  But  it  is  obvious  that  there  is  much 
painstaking  work  to  do  particularly  with  regard  to  pressure 
and  depth.  I  will  come  back  to  this  point  later,  in  connection 
with  modified  ray  theory. 

Computations 

This  subject  has  reached  a  satisfactorily  advanced  stage,  in 
consonance  with  the  enormous  and  rapid  machines  that  are  now 
available . 

"*1  have  a  very  strong  feeling,  however,  that  now  is  the  time  to 
stand  back  and  try  to  fit  the  computational  aspects  of  ray 
tracing  into  a  sensible  perspective.  A  good  example  is  the 
use  of  a  mini-computer:  obviously  essential  in  small 
laboratories,  or  on  board  ship,  where  it  is  an  invaluable 
tool  in  conducting  fruitful  experiments  at  sea;  but  equally 
valuable  in  computer  investigations  where  the  ultimate  in 
accuracy  or  detail  is  not  required  (and  would  indeed  be  a 
hindrance,  especially  when  such  things  as  turn-around  time 
and  availability  are  concerned)  such  as  in  O.R.  work,  tests 
of  ray  tracing  sensitivity,  etc.  But  my  general  point  is  that 
we  should  seek  a  rational  balance  between  such  things  as 
computational  sophistication,  required  accuracv,  accuracy  of 
input  data,  computer  availability,  cost,  flexibility  —  and 
even  decay  time  of  scientists'  enthusiasm. 

To  go  back  to  the  beginning:  we  now  have  a  variety  of  options 
with  which  to  fit  the  sound  velocity  profile  —  linear  segments, 
continuous  gradients,  Epstein  segments,  polynomials,  splines, 
and  so  on. 

This  leads  to  the  fundamental  computational  que  vtion  of  whether 
(a)  to  fit  the  profile  with  segments  for  which  there  is  an 


easily  obtained  solution,  or  (b)  to  use  something  like  a  spline 
fit  and  ask  for  a  differential-equation  solution.  The  linear- 
segment  fit  still  seems  the  most  convenient  choice  for  mini¬ 
computers.  But  on  large  computers  the  choice  between  (a) 
and  (b)  when  one  looks  ahead,  say,  for  the  next  five  years 
seems  much  more  problematical.  Perhaps  approach  (a)  is  more 
suited  to  plotting  rays  and  calculating  travel  times,  whereas 
the  differential  equation  approach  yields  more  reliable 
intensities.  It  may  still  be  too  early  to  decide  the  point, 
but  the  point  nevertheless  needs  attention. 

What  does  the  ocean  think  of  all  this?  Possibly  it  feels  a 
bit  like  a  transvestite  The  sound  velocity  xn  a  step-structure 
ocean  is  a  series  of  linear  segments,  separated  by  discontinuities. 
So,  the  sound  velocity  is  not  necessarily  continuous  in  its 
zeroth,  first,  or  second  —  order  derivatives.  But  more  of 
this  later. 

We  did  not  hear  very  much  about  ray-tracing  when  the  sound 
velocity  has  u  bivariate  profile,  although  we  are  told  that  this 
often  is  the  case  in  reality.  From  a  computational  view  point 
there  is  again  the  choice  between  (a)  seeking  simple  ray 
solutions  in  rectangular  sections  corresponding  to  particularly 
chosen  forms  of  profile  and  (b)  seeking  differential-equation 
solutions  in  regions  where  the  velocity  description  is  made 
by  something  like  bicubic  splines.  However,  we  are  not  aware 
of  any  useful  simple  solutions,  so  the  differential  equation 
approach  looks  more  promising  at  the  present  time. 


Finally,  with  regard  to  ray  tracing  computations,  what  is  the 
general  feeling  about  applying  tests  in  the  manner  of  Moler  and 
Solomon's  use  of  Pederson  and  Gordon's  Epstein  profile  results? 

We  hax  recently  launched  a  modest  venture  at  SACLANTCEN  in  this 
direction,  with  beneficial  results.  Has  anayone's  exprrience 
suggested  using  Epstein  profiles  different  from  that  used  by 
Moler  and  Solomon?  It  seems  that  such  test  programs  would  be 
invaluable  in  the  evaluation  of  ray-tracing  programs,  regarding  not 
only  their  accuracy,  but  also  their  convenience,  speed  and  cost. 
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Experiments 


Due  to  the  unclassified  nature  of  our  discussions  in  the  last  few 
days,  there  may  be  many  more  experiments  than  we  have  heard  about. 
Nevertheless,  I  feel  that  the  experimental  validation  of  ray 
tracing  has  hardly  started.  This  is  a  pity,  since  the  challenge 
of  obtaining  consistent  and  convincing  explanations  of 
experimental  results  is  the  mainspring  of  modern  science  —  it 
would  have  very  beneficial  effects  on  the  development  of  new 
theories  and  computational  techniques.  Of  course,  one  should 
continually  remind  onself  of  the  painstaking  care  required  and, 
with  experiments  at  sea,  the  sheer  physical  difficulty  of  even 
performing  them.  And  overall,  it  seemed  to  me,  "the  agreement 
between  theory  and  experiment  ’'as  encouraging". 


But  for  the  future,  one  can  hope  for  more  and  better  experiments. 
It  would  be  particularly  useful  to  have  many  more  experiments 
of  the  joint  acoustic-oceanographic  type,  such  as  the 
experiment  just  described  by  T.D.  Allan,  in  which  d^nse  spatial 
and  temporal  sampling  of  the  temperature  structure  of  the  ocean 
is  combined  with  acoustic  propagation  experiments,  and  the 
acoustic  propagation  through  known  fronts  just  described 
by  M.J,  Daintith. 

Extensions  to  Ray  Tracing 

There  appears  to  have  been  a  considerable  increase  in 
understanding  of  the  behaviour  of  sound  fields  in  the  region 
of  caustics,  turning  points,  and  shadow  boundaries.  Some 
aspects  of  modified  raj'  analysis  are  positively  sedative. 
However,  these  topics  are  excellent  example^  of  where  the  acid 
test  of  comparison  with  experiments  is  essential  for  future 
progress  to  occur j  experiments  not  only  in  sound  propagation  but 
also  in  the  determination  of  the  sound  velocity  itself  at  depth. 


Theory 

I  started  by  admitting  that  ray  theory  is  an  old  subject,  but 
the  papers  on  the  Riesz  Potential  and  Hamiltonian  methods 
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demonstrated  that  fresh  view  points  can  be  taken  which  may  offer 
the  benefit  of  techniques  adopted  in  other  disciplines.  They 
certainly  merit  very  careful  consideration.  The  demonstration 
of  the  equivalence  of  the  ray  and  wave  approaches  promises 
advances  in  the  description  of  propagation  in  a  medium  with 
statistically  defined  properties. 

It  was  encouraging  to  hear  of  efforts  being  made  to  put 
quantitative  limits  on  the  oft -quoted  conditions  for  the 
validity  of  ray  theory.  Obviously  a  great  deal  more  needs  to 
be  done  in  this  direction,  with  particular  attention  paid  to 
the  obliquity  of  the  ray,  as  well  as  the  scale  size  of  the 
irregularities  in  terms  of  acoustic  wavelength. 

A  related  topic  is  the  validity  of  ray  tracing  at  the 
discontinuous  sheets  separating  the  layers  of  a  step- structure 
ocean.  Of  course,  these  discontinuities  are  physically  large 
gradients,  which  can  be  viewed  as  discontinuities  under 
certain  conditions  of  sheet  thickness  to  wavelength  ratio  and 
ray  angle.  But  it  is  important  to  know  what  these  conditions 
are;  and  what  happens  when  the  conditions  are  violated. 

We  heard  a  little  about  what  we  might  expect  when  a  sound- 
velocity  profile,  or  other  oceanographic  feature,  is  specified 
statistically.  This  is  an  important  avenue  that  has  hardly 
been  entered.  Though  difficult,  statistical  approaches  must 
be  used  to  deal  honestly  with  the  real  ocean.  Perhaps  the  most 
promising  first  step  is  to  examine  the  sensitivity  of  ray  tracing 
to  perturbations  in  the  gradient  or  other  features  of  the 
profile,  and  then  to  extend  this  to  statistically  specified 
profiles,  either  by  analytical  or  numerical  means. 


Applications 

The  application  of  ray  tracing  to  reverberation  modelling  is 
obviously  receiving  considerable  attention,  with  encouraging 
agreement  with  observations. 


It  was  salutary  to  be  reminded  that  the  most  important  channel  of 
our  work  lands  up  on  the  desk  of  the  sonar  systems  designer,  am. 
extremely  useful  to  have  stated  with  such  clarity  his  requirements . 


CONCLUDING  GENERAL  DISCUSSION 


There  was  very  little  reaction  by  the  meeting  to  th?  oceanographers' 
revelations  concerning  the  step-like  structure  of  the  sound  speed, 

A  plea  was  made  for  more  salinity  data  in  these  investigations. 

It  was  pointed  out  that  it  was  necessary  for  acousticians  to 
specify  the  fineness  of  detail  required  in  oceanographic 
measurements . 

Conventional  Ray  Tracing 

The  considerable  fraction  of  the  discussion  devoted  to  the  various 
aspects  of  conventional  ray  tracing  indicated  a  general  concern  with 
a  need  to  refine  present  methods.  Some  of  the  points  raised  in 
this  connection  were; 

Is  it  possible  to  obtain  reliable  answers  while  still  employing  a 
linear-segment  approximation  to  the  sound-speed  profile,  either  by 
taking  a  sufficiently  large  number  of  segments,  o.  oy  smoothing  the 
output  in  some  suitable  manner? 

There  is  a  possible  danger  in  the  use  of  spline  fits  to  sound-speed 
data  that  it  might  sometimes  introduce  artificial  wiggles  in  the 
sound-speed  profile,  perhaps  leading  to  local  gradients  of  sign 
opposite  to  the  actual  gradient,  for  example.  The  opposing  view 
was  put  that  "splines  are  fine",  provided  that  one  works  with 
sufficient  data  points  and  takes  care  with  the  end  conditions. 

An  alternative  approach  to  making  a  rational  fit  to  sound-speed 
data  was  that  the  fit  should  be  contiguous  up  to,  and  including, 
the  second  derivative,  bub  that  the  correct  morphology  should  be 
retained  —  presumably  judged  on  oceanographic  grounds. 

The  possibility  was  mentioned  of  finding  the  eigenrays  by  thp  direct 
numerical  application  of  Fermat's  principle,  as  an  alternative  to 
the  "shooting"  method. 
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A  variety  of  profiles,  for  which  the’-e  are  analytical  solutions 
is  available  from  L.P  Solomon,  Tetr.  Tech-,  Inc., 

1911  Fort  Meyer  Drive,  Suite  601,  Arlington,  Virginia  22209,  U.S 
These  could  be  useful  for  testing  ray  tracing  programs. 

Range-Dependent  Ray  Tracing 

A  recurring  theme  during  the  discussion  was  the  need  to  extend 
present  ray  tracing  capability  to  bivariace  sound-speed  profiles, 
and  eventually  to  trivariate  profiles.  But  there  were  no  clear 
ideas  put  forward  as  to  how  this  could  best  be  achieved. 

Current  methods  include  simply  dividing  the  range  into  blocks, 
each  block  having  a  sound  speed  that  depends  on  depth  only,  and 
then  proceeding  with  the  usual  univariate  methods  within  each 
block.  Another  technique  is  to  divide  the  range  between  given 
profiles  into  triangular  sections,  the  ray  paths  within  each 
section  then  being  circular  arcs.  Both  these  techniques  suffer 
from  implying  oceanographically  unacceptable  sound-speed 
structures  in  regions  between  the  given  profiles. 

Mention  was  made  of  the  Hudson  Laboratories  technique  employing  a 
special  form  of  double  Taylor  series  expansion  (linear  variation 
in  range,  linear  plus  curvature  term  in  depth)  specifying  the 
sound  speed  in  the  region  between  two  given  profiles.  There  is 
also  the  possibility  of  representing  the  bivariate  sound  speed 
in  terms  of  doublj  -  ;ubic  spline  functions.  Ir.  both  these  cases 
the  ray  tracing  could  then  be  accomplished  by  numerically  solving 
the  ray  differential  equations. 

Statistical  Aspects  a.  1  Profile  Sensitivity 

A  difficult  problem,  but  one  requiring  urgent  attention,  is  that 
of  alleging  for  the  effect  of  the  variability  of  the  sound-speed 
profile.  This  can  be  viewed  as  a  problem  of  the  sensitivity  of 
ray  tracing  to  certain  perturbations  o  '  the  so.  /.d-speed  profi  le. 
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Or  it  can  be  viewed  as  a  statistic.!  problem,  requiring  tbo 
statistics  of  the  output  of  ray  tracing  given  the  statistics  of 
the  input  profile(s). 
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The  sensitivity  problem  is  impo'  tu  .t  when  considering  the  adequacy 
of  proposed  methods  of  curve  fitting  to  the  given  sound-speed  data 
points . 

A  major  difticultv  with  the  statistical  problem  is  that  of  posing 
the  problem  properly.  For  example,  if  one  starts  w'ith  a  mean 
profile  (obtained  by  averaging  a  large  collection  of  profiles) 
this  profile  may  be  devoid  of  all  the  features  (such  as  layer 
depth)  which  are  known  to  affect  the  sound  propagation  drastically. 
Therefore,  some  method  should  be  fo'<nd  by  which  the  essential 
character  of  the  profiles  is  preserved.  In  other  words,  one 
should  work  with  a  "typical”  profile,  rather  than  with  a 
strictly  "mean”  profile. 

But  more  significant  than  working  with  mean  values,  the  second-order 
statistics  of  the  rav  tracing  output  are  ai.  important  measure  of  the 
vai lability  caused  by  variations  of  the  sound-speed  profiles  as 
inputs  to  ray  tracing  programs. 

General  Points 

The  opinion  was  expressed  that  insufficient  attention  had  been  paid 
to  the  final  objectives  of  ray  tracing  in  an  ASW  context.  If  we 
are  not  approaching  anything  of  value,  we  might  at  well  stop  now 
and  turn  our  attention  to  potentially  more  fruitful  subjects, 
such  r.s  loudspeaker  design1. 

Two  drastic  alternatives  to  the  present  highly  computer-orientated 
approach  to  the  solution  of  underwater  sound  propagation  problems 
were  proposed.  One  was  to  replace  computers  bj  mathematicians  who 
would  be  cheaper,  and  whose  task  would  hr  .o  develop  al.ternati  ve 
and  more  amenable  theoretical  approaches.  The  other  alternative 
was  to  de-emphasize  computers  and  mathematicians,  and  to 
accentuate  m  compensation  experiments  at  sea.  In  other  words, 
it  might  be  easier  and  more  reliable  to  use  the  Oceanic  analogue 
c.ompu  .  er . 
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There  is,  naturally,  a  bias  of  American  interests  towards  low- 
frequency,  long-range  propagation,  whereas  SACLANTCEN  interest 
is  concentrated  more  on  relatively  short  ranges  at  sonar 
frequencies.  Thus  the  outcome  of  modified  ray  theory  appears 
to  be  of  greater  interest  to  the  former  than  to  the  latter; 
although  some  intriguing  discrepancies  between  experiment  and 
conventional  ray  theory  at  the  shorter  ranges  need  to  be 
scrutinized  in  the  1 » ght  of  modified  ray  theory. 


